Rapid Training of Hamiltonian Graph Networks Using Random Features
Authors: Atamert Rahma
Chinmay Datar
Ana Cukarska
Felix Dietrich
Venue: ICLR 2026 Poster (2026)
Keywords: Graph neural networks
physics-informed machine learning
random feature methods
gradient-descent-free training
Hamiltonian neural network
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Our random feature-based training algorithm for Hamiltonian graph neural networks achieves comparable accuracy while being 100–1000× faster than traditional gradient-descent-based iterative training.

Abstract: Learning dynamical systems that respect physical symmetries and constraints remains a fundamental challenge in data-driven modeling. Integrating physical laws with graph neural networks facilitates principled modeling of complex N-body dynamics and yields accurate and permutation-invariant models. However, training graph neural networks with iterative, gradient-descent-based optimization algorithms (e.g., Adam, RMSProp, LBFGS) often leads to slow training, especially for large, complex systems. In comparison to 15 different optimizers, we demonstrate that Hamiltonian Graph Networks (HGN) can be trained 150-600× faster - but with comparable accuracy - by replacing iterative optimization with random feature-based parameter construction. We show robust performance in diverse simulations, including N-body mass-spring and molecular dynamics systems in up to $3$ dimensions and 10,000 particles with different geometries, while retaining essential physical invariances with respect to permutation, rotation, and translation. Our proposed approach is benchmarked using a NeurIPS 2022 Datasets and Benchmarks Track publication to further demonstrate its versatility. We reveal that even when trained on minimal 8-node systems, the model can generalize in a zero-shot manner to systems as large as 4096 nodes without retraining. Our work challenges the dominance of iterative gradient-descent-based optimization algorithms for training neural network models for physical systems.

Published as a conference paper at ICLR 2026
RAPID TRAINING OF HAMILTONIAN GRAPH NETWORKS USING RANDOM FEATURES
Atamert Rahmal Chinmay Datar Ana Cukarska 12 Felix Dietrich 1.2,4
Tecchncca University of Munich
Munich Center for Machine Learning
TUM Institute for Advanced Study
Munich Data Science Institute
ABSTRACT
Learning dynamical systems that respect physical symmetries and constraints re- mains fundamental challenge in data-driven modeling. Integrating physical laws with graph neural networks facilitates principled modeling of complex N-hody dynamics and yields accurate and permutation-invarianl models. However, training optimization algo-
graph neural networks with iterative, gradient-descen rithms (e.g., Adam, RMSProp. LBFGS) often leads slow training, especially for large, complex systems. In comparison to 15 different optimizers, we demonstrate that Hamiltonian Graph Networks (HGN) can be trained 150-600x faster-but with comparable accuracy-by replacing iterative optimization with random feature- based parameter construction. We show robust performance in diverse simulations, including N-body mass-spring and molecular systems in up to 3 dimensions and 10,000 particles with different geometries, while retaining essential physical invari- ances with respect to permutation, rotation, and translation. Our proposed approach is benchmarked using a NeurIPS 2022 Datasets and Benchmarks Track publication further demonstrate its versatility. We reveal that even when trained on minimal 8-node systems, the model can generalize in a zero-shot manner to systems as large as 4096 nodes without retraining. Our work challenges the dominance of iterative cradient-descent- based optimization algorithms for training neural network models for physical systems.
Kode posttions
Layer trained by random feature sampling
pocttoon (predined)
Φv
V
global pooling
W
H
Layer trained by iqaates optimization
E
go
local pooling
m
G
aocoao
000000 00000000000oooo
Feature processing
Graph processing
Time
Figure 1: We propose efficient training method for Hamiltonian graph networks using random feature sampling and linear solvers (left, also see Figure The HGN captures ground truth dynamics of physical systems (shown: chain of 10 nodes, trained 5) and trains up to 600x faster than State-OffThe-Art (SOTA) optimizers.
INTRODUCTION
Learning from data requires careful design in several key areas: the data, the model, and the training processes. To enable the model to generalize beyond the training set, it is important to incorporate set of inductive biases into these processes. When approximating physical systems, is beneficial to include physical priors to accurately capture the system's 's characteristics, including its dynamics and the fundamental physical laws (Tenenbaum et al.) 20000; Chang et al.) 2016 Watters et al. 2017; de Avila Belbute- Peres et al.) 2
Sharma & Fink (2225). Consequently, many architectural designs are rooted modeling physical frameworks, such as Hamiltonian mechanics (Bertalan et al.
2019 Greydanus al.) 2019), Lagrangian mechanics (Cranmer et al. 2019 Lutter et al. 2019
"Corresponding author, atamert.rahma tum.de

---

Published as a conference paper at ICLR 2026
-ber-Bloebaum& & Offen 2023). Nort-Hamiltonian systems (Desai et al. 2021 Roth et al. 2025).
and GENERIC (Heráández et al. 2021 Lee et al. 2021 Zhang et al. 2022b; Gruber et al. 2025
Graph networks have useful inductive biases such as locality and permutation invariance (ooso
ctt al. 2024p. which are desirable for many interconnected, complex systems observed in nature.
Thus, for many applications in natural sciences, graph network model is suitable choice. The
key aspects of modern graph networks include neural message passing (Gilmer et al.] (017) and encoding additional local information into the system (Corso et al.) 2024, Schlichtkrull al. 2018
Brookschmidt, 2020). In physics, graph networks have been employed to analyze data from the Large Hadron Collider (DeZoort et al. 2023). model mechanical systems (Zhao et al. 2024). and fluid
dynamics (Xue et al.j Peng et al. 2023 et al. (0024).
Efficient and robust training of graph networks on large systems for natural and life sciences is an active area of research. Despite the advantages of using graphs for physical N-body systems, their training reportedly slow due pradient-descent-baseee iterative optimization (Kose Shen
2023 Shukla et al., 2022; Vignac et al.) [020 Kumar et al. 2023 Zhao & Cheah, 2025 Marino et al. 2025). These challenges become even more pronounced when a numerical integrator is incorporated
into the model architecture (Xiong et al. 202T). Furthermore, physics-informed models are often
sensitive to the selection of hyperparameters (Shukla et al. 2022), which increases the challenges
posed by slow iterative training.
Recently random feature (RF) networks have been shown to be promising for approximating
physical systems (Fabiani et al.
Datar et al.) 2024 Rahma et al., 2024 Bolager et al.. 2024 However, to the best of our knowledge, random features have been used train graph networks for modeling physical systems. Recent work RF-Hamiltonian
2021
Fabiani et al. 2025 Galaris et al.
2022
neural networks (RF-HNNs) is promising (Rahma et al. 2024), where RF-HNNs are trained without using iterative algorithms. The authors demonstrate very low approximation errors, but only for very small systems and without leveraging the graph structure. In this paper, we introduce an efficient and accurate training method that utilizes random features for Hamiltonian Graph Networks (RF-HGNs. see Figure T). Our main contributions are as follows.
We introduce Random Feature Hamiltonian Graph Networks, combining random sam- pling with graph-based physics-informed models for the first time, and show how one can incorporate translation, rotation, and permutation invariance as well as knowledge about the physical system (see Section 3
We provide much faster and more efficient alternative to gradient-descent-based iterative optimization algorithms for training that avoids challenges related slow, non- convex optimization and vanishing or exploding gradients (see Section
We perform comprehensive optimizer comparison with 15 different optimizers and demonstrate that random feature-based parameter construction offers up to 600 times faster training for HGNs, without sacrificing predictive performance (see Section [.1). The demonstrations are performed on challenging benchmark problems, including mass-spring and molecular dynamics systems.
We demonstrate strong zero-shot generalization, with models trained on graphs with tens of nodes accurately predicting dynamics on graphs with thousands of nodes (see Section 4.2).
2 RELATED WORK
Training graph networks: graph structure allows for modeling a wide range of processes due to structural properties or underlying causal relationships. For many problems, the best way to achieve good performance is by using a large, high-quality dataset to train large model. In such settings, training can be significantly slowed down due to the computational effort needed for backpropagation. Improvements can be achieved with specific sampling methods for the training data (Nagarajan &
Raghunathan| 2023 Zhang et al. 2022a Zhou et al. 2022 Kaler et al. 2022
Zhang et al. 202 1
Lin et al. (2200). graph coarsening (Hashemi et al.] 2024| Kumar et al. 2023 Jin et al. 2022a b
Bravo Hermsdorff & Gunderson
2019), or hardware acceleration Shao et al. 2024 Gupta et al.
2024 Zhu et al. 2025 Wan et al. 2023 Yang et al.) 2022 Kaler et al.) [022] Wang et al. 2021) Cai
al.
2021 Lin et al. 2020). Nevertheless, graph networks for physics still face challenges during
2

---

Published as a conference paper at ICLR 2026
training due to need for high accuracy in the dataset, irregular memory access, load imbalance during backpropagation (Shukla et al. (022), and hyperparameter tuning (Schmidt et al. (021).
Graph networks for physies: A notable advantage of graph-baser models is that they are tractable for nigh-dimensional data, assuming that the graph connectivity remains sparse, which is suitable for many physical systems learned in a data-driven way. Recent work has incorporated graph neural
networks into their model architectures for approximating physical systems (Pfaff et al.) 202T]
Gonzalz- et al.) 2019; Sanchez-Goozalez et al.] 2020: Tierz et al. 2025] Varghese et al.
2025 Bhattoo al., 2022, Thangamuthu al., 2022). However, training graph network for Very large number of nodes is challenging; one possible remedy is to partition a large graph and enable information exchange between partitions, which are trained individually (Nabian et al. (0024). Approaches for graph networks for Hamiltonian and Lagrangian systems are used by Thangamuthu
et al. p: Bhattoo et al. 2022); Bishnoi et al. (2023); these models are typically trained with the Adam optimizer and applied to N-body systems. Other work addresses the issue of long-range
information loss in large graphs by adding physics- based connections (Yu et al., 2025), yielding
better predictions but not addressing training difficulties that might arise.
Random features for graph networks: Random features originated with the idea of using perceptron by RRoeenlaatt (1962) and gained traction after theoretical contributions established that
they can lead to accurate approximations (Johnson & Lindenstrauss 1984] [arron 1993 Rahimi
& Recht 2007) at low computational cost. As the machine learning community gains a better
understanding of random features (Rahimi & Recht (2008; Bolager et al. 2023 Fabiani 2024). many
new variants are being explored Zozoulenko et al. [022 |Bolager et al. 2022]] Datar et al. 2024
Rahma et al., 2024). An innovative approach for graph classification problems used random features
approach and demonstrated competitive accuracy on large classification datasets with training time of only few seconds minutes (Gallicchio & Micheli 2020). Such an approach is related to
echo state graph networks (Gallicchio & Micheli] [010: Wang et al., (023). Recent work has also developed graph random features enabling kernel methods on large graphs (Choromanskij [2023] Reid
iet al. 2023b a.. leading to a notable reduction of the cubic time-complexity for kernel learning.
3
METHOD
Problem setup: In this study, we aim to efficiently learn the Hamiltonian of dynamical system from observed phase space trajectories, while exploiting the underlying graph structure and incorporating relevant physical invariances into the model. We consider a target Hamiltonian for an N-body system on R2d-N the Fuelidean phase-space of dimension 2d N E N, where d is the spatial dimension. We denote the generalized position and momentum vectors by q.p E R with q, pl Rd denoting the 7th particle's ' state. We denote by ă the time derivatives of a trajectory z(t): : IR + & for & e N. The Hamiltonian is valued function H 2d - R that describes the system dynamics in
the phase-space through Hamilton's 's equations (Hamilton 1834 1835) given by
q P
= JVH(g.P), J=
0 I -I
E ER (g(2d-N)x(2d-1))
0
where / [gdd-Nx(ddNN) is the identity matrix. We summarize the notation in Appendix A
We focus on N-body systems in this work, where the graph representation is naturally available, c.g., chain of masses connected via springs. Given dy, dg N, we write the system with N nodes as a graph G = (V, E') with node feature set V = {o; E dd ii=l == L..... ......N] N} and an edge feature set F. (ci E Rd Vi,j such that AGj = 1}. where A NNN is the symmetric adjacency matrix that encodes the node connectivity information. We parametrize the Hamiltonian function with a graph neural network, and then use the trained network to simulate the physical system by integrating Equation (I) with the symplectic Störmer- Verlet integrator ((Hairer et al. 2003), also see (Offen & Ober- Bloebaul ). In contrast to previous work on Hamiltonian Neural Networks (Bertalan et al. 2019 Greydanus et al. 2019: Dierkes et al. (023). train our networks through random feature sampling algorithms rather than iterative, gradient based optimization. Thus, call our approach "rradient-descent-freee.
Graph representation:
3

---

Published as a conference paper at ICLR 2026
3.1 ENCODING INVARIANCES
The systems we consider are translation-, permutation-, and rotation-invariant i.e., when the whole system is shifted, permuted, or rotated, the Hamiltonian stays constant. With the general Hamiltonian formulation, however, cannot guarantee such invariances by design (Du et al. 2022). To construct such invariant representations, we introduce transformed coordinates q. P E Re derived from the original phase-space coordinates q. IRdNN defined in arbitrary reference frame.
Translation-invariant representation: To make the position representation (and consequently the system representation) translation invariant, normalize the positions by subtracting the mean di - 4: t Σ q.. We do not make the generalized momenta p translation invariant, as shifting the momenta would change the total energy of the system in N-body systems, for instance, when the kinetic energy depends on the norm of p.
Bermutatinniivaaiinne The graph structure and appropriate message passing algorithms inherently provide us with a system representation that is invariant with respect to node index permutation.
Rotation-invariant representation: Starting from the translationiinvarii representation, we then perform an- other transformation make the final representation also rotation-ing potation-nvvrrian. invariant. Here, we explain how to encode a rotation-invariant representation for single-body system N and spatial dimension for brevity.
10-3
Rotation
Train
Test
invariance
MSE
10
10-11
Translation-
100
invariance
11
Using true it
10.9
RF-HGG Boon-nuar
10.8
RF HF-HEN
We construct a new representation in local orthonormal basis starting from the original coordinates p, q E R² To construct the basis, we choose the first basis vector as 61 [Ivill E R.?
Zero-shot generalization
10.7
10.6
10.5
10.4
10.3
Time step
where ||-Il denotes the 12 norm. We construct a second basis vector Rei E RR2, which orthonormal to Q1 IR² and obtained by simply rotat- ing 4i by 90º using rotation matrix R. We then define orthonormal matrix with the two basis vectors as
Figure 2: Illustration train and test N-body system posi- tions showcasing the RF-HGN's translation- and rotation- invariance, and its zero-shot generalization capability, vali- dated by conserved Hamiltonian and low trajectory predic- tion errors for the test data (see Appendix B. for details).
B B=
62). Finally, the rotation-invariant coordinates for any body (L....,,N . are q: B3 g:-
Given fixed first point, our procedure yields a rotation-invariant representation. One can uniquely identify the first point, independent of node ordering or orientation, as the one closest to the mean q. In case of ties, we select the point with the smallest angle relative to the first coordinate axis centered q. The same procedure is applied to obtain iaaaiioninvvaiann representations of the momenta. In higher dimensions, when d > 2 and N 1, we follow similar procedure, but construct the orthogonal bases using Sram-Sch orthogonalizatton instead (see Appendix [B.2]. Figure 2 demonstrates how translating the N-body system, rotating it, and even adding new nodes without re-training (erro-shot-eenerrlizatoon,) yields low trajectory errors while conserving the Hamiltonian.
3.2 MODEL
We now describe each component of the architecture of the Hamiltonian Graph Network (HGN) in detail (see Figure B). Please refer to Appendix B.3 for details on the forward pass.
3.2.1 NODE AND EDGE ENCODING
Node features: For an N-body system with translation- and rotation-invarian position and momenta representations, we define node features as v; l I IRdv, where dy 2d, for each i E (1,..., NJ. We define the set V = {ra i 1,..., N) that collects all node encodings.
4

---

Published as a conference paper at ICLR 2026
Φ
h
ΦM
h V
global pooling
W
Hl
h
local
m
Ø 25
pooling
Figure 3: Hundom-feature Hamiltonian graph neural network architecture. Left (green box): Con- struction of node and edge encodings have and h from translation and rotation invariant position and momenta p representations of an N-body system. Right (orange box): Construction of a global encoding for the graph using message passing. In RF-HGN, dense layers (blue) are constructed with random features, and linear layer weights (red) are optimized by solving a linear problem.
Edge features: We define the edge features for each edge (i.j) with [( --);;;; ; 1b - -glll'' R&s, where dg = d + 1 and take the absolute value ((q: - 43) ) the molecular system examples. We collect all the edge feature encodings in the set E' == {eis j and i,j E [N.....N].. In order to reduce the memory and computation costs, we define canonical direction by always computing edge features from higher to lower-indexed nodes, such that each edge is represented only once and set (ejs = egy)- We use the relative displacement vector Oj and its norm to represent the direction and distance between connected nodes in the system, order to capture local geometric structure and pairwise interaction properties.
Cjj
in
Input encoding: The constructed node and edge features are then encoded via separate dense layers,
NY (v:) oWW + bv) E Rd
Ye; E V, and
(2)
h1 )n(e))) o Ggj bε) E Rd
Vei E F,
where Φv RR - Rd and dE dR - Rddh are outputs dense layers that encode the node and edge features, respectively, with weights Wy Rd xdy WW Rd xdæ and biases by by R. We denote the activation function (here, softplus or gelu) by a. Using the symmetric edge features described earlier avoids duplicate memory and computation overhead in the input encoding as well, since we only compute the encoding h of each undirected edge feature esj where > and use the same encoding for both directions h and E
3.2.2 MESSAGE PASSING AND FINAL LAYER
We perform bidirectional message passing along edges (i. where A.j 1, allowing nodes to aggregate information from their local neighborhoods.
Message construction: Messages are constructed from the encodings of source and edge nodes via dense layer ΦM 2d - Rd as
h
ΦM
WM
h
bag
E gd
with weights WM dd xR and biases bM RdM,
Message passing (local pooling): Each node aggregates incoming messages using permutation- invariant operation (here, summation) m, EieN, hy where N; is the set of neighbors of node
(source of incoming edges j where Ay 1).
Graph-level representation (global pooling): All node emheddings and aggregated messages are pooled to form a global encoding the network, such that
Thy mj
hG Σ
E Rd,, where d2 = d, dss-
(5)
Linear layer: The graph representation linearly mapped to scalar value that approximates the conserved value (energy) of the system, such that = W,, bG + b,. where W,, e gd and by, E R denote weights and bias of the linear layer, respectively. Without loss of generality, omit the bias
5

---

Published as a conference paper at ICLR 2026
term, it acts only as integration constant and does not affect the dynamics oM and ØH og öp which our primary interest. We assume this constant is known to align the model's conserved quantity with the true Hamiltonian H in all examples.
3.3 TRAINING
We now describe the central idea of this paper our training algorithm, where we use random feature sampling techniques instead of gradient-descent-baseee iterative optimization algorithms. In particular, we discuss how to compute dense and linear layer parameters of the network.
Dense layer parameters: We compute the weights and biases Wy, We, WM. by bg. bus of all the dense layers (øv, Φg, ΦMr) using random sampling algorithms. Specifically, we use two sampling approaches here: Extreme Learning Machines (ELM) (Schmidt et al. 1992 Pao & Takefuji 1992, Huang et al.] 2004 2006 Rahimi & Recht 2008; ZZang et al. 2012 Leung et al. 20199) and the "SSample Where Matters" (SWIM) algorithm (see (Bolager et 2023)) for unsupervised learning problems (see (Rahma al. 2024 Datar et al. 2024)). As an illustrative example, describe how to compute the parameters of the dense layer øy and use the notation from Equation ( The parameters of the other dense layers sampled analogously.
The (ELM) RF-HGN approach is data-agnostic. The weights Wy are sampled from the standard normal distribution, and biases by from the standard uniform distribution. The (SWIM) RF-HGN ap- proach is data-driven. The network parameters are computed from pairs selected uniformly at random from input data points Fú, where the point coordinates correspond to dense layer inputs. The weight and the bias of the hth neuron in the dense layer are constructed using the input data pair (x)
chosen uniformly at random from all possible pairs, so that w; = $1 (.x (2)
x2 (2) ) [z]
(1)
and b, - (w,, a (1) $2. Here, (w,,b.) are the weight and bias of the hh neuron, and ($1. $2 are
constants depending on the activation function used in the dense layer (see (Bolager et al. 2023 for an analysis).
Linear layer parameters: After sampling all the dense layer parameters, we compute the optimal parameters for the linear output layer of the network by computing the least squares solution (see (Rahma et al. 2024), but also related work (Bertalan et al., (019)). For an N-body system, we denote a single input to the RF-HGN by y E Rzav and the output of the global pooling layer by
((v) E R2 and the total number of input data points by M. The linear system that approximately satisfies Hamilton's 's equations is then
(
(p(UM)
T (2)))
L bz
1 J
UM H(yo)]
(6)
WEBEd-N.MM-
N (1)
Equation solved for the linear layer parameters Wr, and by using 12 regularization. The regularization constant in most examples was chosen very small (see Appendix D). We assume the true Hamiltonian value H(yo) to be known for a single data point to fix the integration constant br. We assume there is external force acting on the system during training, such that the total energy is conserved. However, we can easily add an external force while evaluating the trajectory during inference. In the computational experiments, given time derivatives æ. We demonstrate training the model purely from time series data as part of the benchmark experiments (Section 4.4).. Equation (6) results in a convex optimization problem, [WT_)bL+1] b arg ming, [|ZA] u/I which can be solved using efficient least-squares
mostly train with explicitly
algorithms (Meng et al. 2014b).
Runtime and memory complexity: During training, the sampling of dense layer parameters and gradient computation are fast, with the primary run-time bottleneck being the least squares solve (Equation (6). Assuming d, & K = 2d N . (which is always the case in our experiments), the total run-time complexity is OKdr,)). This highlights an important feature of our approach: training time scales linearly with data size M, the number of particles N, and the spatial dimension d, given fixed settings for other variables. The memory complexity during training O ((MN. and thus also scales linearly with the number of edges, and dataset size (see (Bolager et al. 2023 and Appendix B.4 for details).

---

Published as a conference paper at ICLR 2026
4 COMPUTATIONAL EXPERIMENTS
We evaluate our method on mass-spring and molecular dynamics systems with two and three degrees of freedom (2D and 3D), as illustrated in Figure]4 We provide further details the used datasets, setup, and hardware in Appendix|C [C) Appendix [ and Appendix [ respectively. We additionally discuss hyperparameter tuning and ablation studies of increasing feature widths and number of message passes in Appendix F test robustness against noise in Appendix [H.4] demonstrate batch-
wise training Appendix H.S and provide further random feature benchmarks Appendix H.6 For
Tuue: 2prd E R for m E N. we define the relative eπror as Tuue T 2// [
4.1 BENCHMARKING AGAINST SOTA OPTIMIZERS
Table I: Results of training the HGN architecture for the 3D lattice system (see Figure 4 (a)) with different optimizers, Results show mean (min, max) over three runs on the same GPU hardware.
Optimizer
Test MSE
Train time [s]
Speed-up
RE-HGN (ours)
8.95e-5 (6. (6.96e-5, 1.13c-4)
0.16 (0.13, 0.22)
LBFGS (Liu Nocedal [989]
3.56e:5 (1.21e-5, 5)
(23.711 23.95)
148.96x
Rprop (Rrcumnter & Bralin 1993
9.59e (7.49e-5, 2.63e-3)
30.84 (30.74, 30.94)
192.62x
RMSpepp [Telemam & Tunton 2012
1.9e (2.55e 3)
91.62 (91.13, 92.42)
S72.24x
Adam (Kingma
2015
2.90s (4. 11e-5 3.53e.3)
91.64 (89.97, 92.66)
572.37x
SGD SGD+monentum Shiikkevrr ot al.
4.23e-3 (3.81e-3, 4.06i-3)
91.65 (91.14, 92.20)
572.43x
2013
SGD Rubbins [951]
2.36e-2 (1.85e-2, 2.86e-2)
91.75 (91.51, 91.91)
$73.07x
Adagrad elchhellall
2011
2.58e (2.88e-3, 7.05e 2)
92.03 (91 18,58,88000100 92.47)
574,8 84x
Adam
Hilfer
2019
2.91e:3 (4.30e-5, 8,53e-3)
92. (91 86, 92.31)
S75.59x
Adamux
Kingia & 2015
1.5e 31.5ee 4.332e
92.33 (92.07, 92.69)
$76.68x
Adadelta
/Zaler 2012)
(1.49e3, 1.96e-2)
92.60 (92.25 93.07)
578.39x
Radam Liu cr 2021
1.69-- (5.36e-5, 4.75c-3)
93.00 (92.75, 93.42)
580.88x
Dozal 2015
9.1ee4 (4.00ee,,,
93.42 (92.88, 93.80)
583.54x
Averged
SGD
Lower al
2019
2.36e 21.85e 2,86e 2)
94.50 (94.01, 94.78)
$90,26x
Adafactor (Shazeer & Ster
2.1s:3 (1.06e3,) 4.71e.3)
96.36 (95.67, 96.88)
601.92x
The goal of this experiment is to demonstrate the efficiency of our train- ing approach in comparison with the conventional training methods that rely on SOTA iterative optimization algorithms. To this end, we consider all of the existing optimizers available in Py Torch (Paszke et al.) 2219)) as the current SOTA iterative training proce- dures. In this experiment, the target function is the Hamiltonian of gen- eralized N₂ x Ny-body lattice mass- spring system with a spatial dimen- sion d = 3 given by Equation (C.8). Table lists the results of training the HGN architecture for the lattice sys- tem with different optimizers, sorted by training time. The hyperparameters are tuned for each optimizer separately, and early stopping was used for all iterative approaches. Our proposed training method significantly outperforms all iterative approaches in terms of training time by factor of 148 up to 601, and is only slightly less accurate compared to the LBFGS method (a second-order optimizer). Section 4.4 includes comparisons for other graph network architectures on benchmark datasets, where we observed similar results.
(a)
(b)
(c)
Figure 4: Graphs considered in the experiments: (a) 3D lattice (nodes arranged on a 2D grid, moving in 3D space- see Section 4.1 and Section 4.2). (b) an open chain (nodes moving in 2D space interacting through Lennard-lones potential (nodes moving in 2D space with dynamic edges see section Section (.22 and (d) 2D closed chain (nodes moving in 2D space see Section 4.4
Section 4.2 and 4.3) (c) molecules
4.2 ZERO-SHOT GENERALIZATION AND COMPARISON OF RANDOM FEATURE METHODS
We now study zero-shot generalization, where we train an RF-HGN on small systems size 2x2, 3x3, and 4x4, and test systems going from 2x2 up to 100x 100. approximate Hamiltonian for much smaller systems with 3x3 nodes and reliably predict it with extremely large systems of size 100x 100 without retraining. A 2x2 system is an edge case where all
Figure 5 shows that
accurately
7

---

Published as a conference paper at ICLR 2026
nodes have only two edges, lacking the nodes with four edges, in the test data, explaining the poor zero-shot generalization.
10
Training Systems
Systems
erfor
2x2
3x3
10x10
50x50
100x100
10
Train system: 2x2
Train system: 3x3
10-1
Train system: 4x4
4x4
10
242
3x3
4x4
5x5
10:10
20x20
50x50
100x100
system
Figure 5: Illustration of accurate zero-shot generalization for 3D lattice (see Figure
4 (a)): Training on
smaller systems (left) enables accurate predictions (right) on extremely large test systems (middle).
using true # ....
(ELM) RF-HIN
(SWIN) RF-HNN
(ELM) RF-HON
(SWIM) RF-HGN
0.2-
02
30-
1014-
9: of node
0.1
9: of node 4
0.0+
10.1
-0.1-
-0.--
20-0-=
10
-0.3
-044
10"10
2000
4000
2000
4002
2000
4000
2002
4000
Time step
Time step
Time step
Time
633
0. node 2048
node 2048
122
10-
MSE
10"
0.0
-0.05
-011
-0.00
435x162
10:13
-0.15
4.1
2000
Time step
Time step
Time step
Time stop
Figure 6: Illustration of position trajectories (first two columns), their true Hamiltonian values (third column), and MSE (fourth column) over time from models trained system with eight nodes on
the 2D open chain (see Figure 4 (b)). Top row: Results from RF-HNN and RF-HGN architectures are visualized along with the ground truth. A system of the same size for training and testing is used (23 nodes). Bottom row: Results from RF-HGN architectures and ground truth for zero-shot test case with system size of 22 nodes (trained on 23 nodes).
We consider another example of an N-body chain mass-spring system in 2D space (eq. C.99 Figure 16 11.166shor shows that by training on a much smaller system with 2 nodes, RF-HGN trained with ELM and SWIM results in very low errors for systems large 22 nodes. This demonstrates robust and strong zero-shot generalization in graph-based architectures, with training using SWIM outperforming ELM by approximately an order of magnitude. Since HNNs are not graph-based architectures, they have to be re-trained for each system, and even after re-training, they perform much poorly in comparison with graph-based architectures (by 1-2 orders of magnitude). We also observe that the graph- based architectures (RF-HGN) are slower to train compared to their respective counterparts trained with fully connected networks (RF-HNN), but are more accurate by 1-2 orders of magnitude, especially for large systems, even without re-training.
Figure 6 (top row) shows that the trajectories evaluated with the RF-HGN accurately match the true ones closely while approximately conserving the Hamiltonian, unlike RF-HNN, The RF-HGN trained with SWIM, in particular, outperforms the one trained with ELM by roughly two orders of magnitude. In the bottom networks. Our (SWIM) RF-HGN robust and again exhibits a low error, while the ELM-trained model exhibits deviations from the true Hamiltonian.
of Figure 6 show zero-shot results, which limited to only the graph
We now consider another Hamiltonian, using the Gennard-Jonee potential (eq. C.10 to investigate generalization properties to different geometries. First, small system with 9 particles was trained
8

---

Published as a conference paper at ICLR 2026
with Adam and the random feature methods ELM and SWIM. The accuracy with ELM was poor (see Figure|HL. H8) and thus the results are omitted in Figure II19 We speculate that the main reason for poor performance is the non-isotropic input variables, which makes normal distributions a bad choice for the weights. The results are given in Table Figure 19 TableHH.32 Table H..33 and Figure [ To evaluate RF-HGN in a more complex scenario, employed dynamic edge indices with cutoff of 2.0, trained our model with 36 particles, and tested with 64 particles to test zero-shot generalization. We visualize the rollout trajectories in Figure H.20 and Figure H.2T and observe that SWIM sampling clearly outperformed ELM, while maintaining slightly worse approximation than the Adam-trained HGN. We note that none of the trainers could reach low approximation errors (~ 10% relative error). To the best of our knowledge, our RF-HGN is the first random feature-based physics- graph network, and can be trained approximately 100 times faster than with the Adam optimizer at comparable accuracy.
Table 2: Molecular dynamics evaluation with 9 particles. Mean squared error (MSH) and relative 1 error (rel. (2) are reported together with the true Hamiltonian over the ground-truth trajectory and the (SWIM) RF-HGN predicted quantity over the rolled-out trajectory.
T=1
T=25000
T=50000
T=74999
T-99999
MSE
1.140e13 13
3.239e-03 03
1.998e 02
4.932 02
8,301e 02
q rel. 12
2.346e 07
3.978e-02
9.717e-02
1.545e-01
2.032e
True
-1.233e401
-1.233c=01
-1.233c+01
-1,233ex01
-1.233es01
Model i
-1.223e401
-1,223e+01
-1,248e+01
-1.214e+01
-1,217e401
4.3 BENCHMARKING WITH REAL-WORLL POTENTIALS WITH INCREASING COMPLEXITY
To further evaluate the applicability of RF-HGN, we experiment with additional potentials from quantum mechanics (anharmonic oscillator Atkins et al. (023) and molecular dynamics (the Morse potential Morse (922)), and list testing results with an unseen initial condition in Table 3 In addition
to using more complicated potential, also now apply an external (gravitational) force node-wise
H. 1. for learned during integration and simulate for a long-time horizon (see Figure 1. 1.
Hamiltonian plots over the trajectory: see Figure H.IT H.13 H.15 snapshot visualization of the predicted trajectories). The results show that more challengigg potentials mnn-llneea forces) than the standard mass-spring (linear force) can also be approximated with the RF-HGN model with reasonable accuracy, compared to the Adam optimizer, still achieving 200-300x speed-ups, even without GPU acceleration. See Table 3 for details, and note that the RE-HGN results are achieved without extensive hyperparameter tuning (as opposed to Adam).
Table 3: Zero-shot (trained with 5 nodes, tested with 10 nodes) test evaluation by solving an unseen initial condition of a 2D chain (see Figure using Hamiltonian graph models trained differently. The system is solved for 10,000 time steps with time step size At = 10 with gravitational force applied during integration, and the last position MSE is reported against the reference solution using the true Hamiltonian. The results are listed for the standard spring, anharmonic oscillator, and the Morse potential in that order.
Posential
(Adam) HGN
(ELM) RF-HGN
(SWIM) RF-HGN
V(r) 10,2 V(r) [0,2 + 1m0
3,875e-03
2.331c-03
3.408c-05
4.562e-02
4.324e-02
5.232e-04
V(r) D(1 explp 2
8.893e-02
7.398e 04
1.218e-03
4.4 BENCHMARKING AGAINST SOTA ARCHITECTURES
The goal of this benehmark is to compare the results of our model with the existing state-of-the-ee graph-based network architectures used to model physical systems. The Adam optimizer (Kingma & Ba Ba||2015) widely regarded as the SOTA optimizer for physics-informe
GNNs (Kumar et al. [023] Thangamuthu al. 2022). is used as the default in our comparisons. We use the dataset and code from Thangamuthu et al. (2022), introduced at the NeurIPS 2022 Datasets and Benchmarks Track. In observe excellent performance the benchmark spring systems (Figure 4 d) with orders of magnitude faster training times while maintaining a comparable accuracy.
line with the other experiments,
Training times are reported in Table In our evaluations, certain specialized architectures such as
9

---

Published as a conference paper at ICLR 2026
Lagrangian Graph Networks (LGN) (Bhattoo al. 2022) occasionally exhibit instability and diverge on test trajectories, whereas our model maintains robust performance. For a detailed problem setup, accuracy comparison, and the architectures, please Appendix G
Table 4: Comparison of training times (in seconds) for RE-HGN optimized using SWIM to existing inysics-informed graph models optimized with Adam on a benchmark dataset from Thangamuthu al. (2022) on 2D closed chain (see Figure (d)).
System size
(SWIM)
RE-
FGNN
FUNODE
GNODE
LGN
LGNN
HGN
HGNN
HGN
406
380.35
2367.37
12534.81
7225.
1288.08
N
3.87
475.24
430.32
2499.04
20536.78 78
6259.58
1370.14
4021.59
5.42
536.27
520.54
2600.31
53148.24
8774.59
1676.78
4380.46
CONCLUSION
We propose a training algorithm for Hamiltonian graph networks via rapid random feature sampling and linear system solvers. Our approach completely avoids slow, iterative pradientdeescen--eee based optimization, which is especially challenging in the graph network and the physics-informet settings. We demonstrate our approach chain, lattice, and molecular systems in up to three spatial dimen- sions, encompassing N-body systems. By incorporating translation, rotation, and index-penn permutation invariances, extend random feature methods to graph-based Hamiltonian network architectures. Compared to 15 optimizer baselines, our method offers dramatic speedups (100x to 1000x) while achieving competitive accuracy in 3D physical systems. Remarkably, training on 3 x 3 systems suffices to accurately predict dynamics in systems of size 100 x 100, demonstrating strong zero-shot generalization capabilities. With this generalization from such small-scale training systems, one can deploy models without needing to re-train on full-scale data, enabling fast prototyping.
Limitations and future work: For very small graphs, the HNN architecture is often faster to train than a graph- based approach, making it better choice than HGN in these cases. Generalization capabilities of the HGN models are typically limited to the same type of graphs, i.c., models trained on chains (edge degrees up to two) cannot be used to predict dynamics of lattices (edge degrees up to four) (Corso et al. 2024). We observed similar challenges when using dynamic edges in the molecular dynamies examples for all optimizers. Therefore, for those examples, one could either use more data or less physically complex model (c.g., an SE(3) graph network that provides SE(3) equivariance by construction (Du et al. (022g) to reduce data requirements. Our approach does not easily generalize to other graph neural network architectures (e.g., convolution or self-attention on the individual node features). future work, we intend to extend this work by employing multiple message passing layers and deeper architectures where random feature boosting might help (Zozoulenko et al. 2025
Ethics statement: We demonstrate that data-driven construction of random features can signifi- cantly outperform many SOTA optimizers in terms of accuracy and training time. The tremendously increased training speeds we report may also speed up the development of nefarious and even danger- ous applications. Similar all HGN and HNN models, for this purpose. However, specific bad intent as well as significant further development would be required for this to happen. We hope that, instead, our work has a profound positive societal impact in the future, because training such accurate models from data is important in many sciences as well as in engincering but has been slow up now, due to the difficulties in training.
specific training method is not designed
Reproducibility statement: We provide details on the used datasets, model setup with hyperpa- rameters, and hardware details in Appendix C. Appendix D and Appendix B respectively for all the numerical experiments discussed in the main text and in the appendix. In the supplementary materials we provide further instructions in a README. mc on how to reproduce all results in the form of tables and plots. Our codebase is publicly available at
ht ht Sttps://gge :ttt::///tiii http:/taa...cceeddeeeeaaa/c///ca
10

---

Published as a conference paper at ICLR 2026
Acknowledgments: We grateful for discussions with Samuel James Newcome, Manish Kumar Mishra, Markus Mühlhäußer, Jonas Schuhmacher, Iryna Burak, Nadiia Derevianko, Qing Sun, and Erik Lien Bolager. A.R. is supported by the BMBF (project AutoMD-AI, A.Ć. is supported by the TUM Georg Nemetschek Institute Artificial Intelligence for the Built World., and C.D. and F.D. supported by the DFG (project no. 468830823), and acknowledge association with DFG-SPP-22. C.D. is partially funded by the Institute for Advanced Study (IAS) the Technical University of Munich. The authors gratefully acknowledge the computational and data resources as well as the support provided by the Leibniz Supercomputing Centre (www.Irz.de).
REFERENCES
Peter William Atkins, R George Ratcliffe, Julio De Paula, and Mark Wormald. Physical chemistry for the life sciences. Oxford University Press, 2023.
A.R. Barron. Universal approximation bounds for superpositions of sigmoidal function. IEEE Transactions on Information Theory, 39(3):3309945, May 1993. ISSN 0018-9448, 1557-9654. doi: 10.1 1..11091/8.25500
Tom Bertalan, Felix Dietrich, Igor Mezic, and Ioannis G. Kevrekidis. On learning hamiltonian systems from data. Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(12). 2019.
Ravinder Bhattoo, Sayan Ranu, and N M Anoop Krishnan. Learning articulated rigid body dynamics with lagrangian graph neural network. In S. Koyejo, S. Mohhame,, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 29789-29800. Curran Associates, Inc., 2022.
Suresh Bishnoi, Ravinder Bhattoo, Jayadeva, Sayan Ranu, and N. M. Anoop Krishnan. Discovering Symbolic Laws Directly from Trajectories with Hamiltonian Graph Neural Networks, July 2023.
Erik L Bolager, Iryna Burak, Chinmay Datar, Qing Sun, and Felix Dietrich. Sampling weights of deep neural networks. In Advances in Neural Information Processing Systems, volume 36, pp. 63075-63116 Curran Associates, Inc., 2023.
Erik Lien Bolager, Ana Cukarska, Iryna Burak, Zahra Monfared, and Felix Dietrich. Gradient- free training of recurrent neural networks. arXiv preprint arXiv:2410.23467, 2024.
Gecia Bravo Hermsdorff and Lee Gunderson. A unifying framework for opectrum-preservvinn graph sparsification and coarsening. Advances in Neural Information Processing Systems, 32, 2019.
Marc Brockschmidt. GNN-FiLM: Graph neural networks with feature-wise linear modulation. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 1144-1152. PMLR, July 2020.
Zhenkun Cai, Xiao Yan, Yidi Wu, Kaihao Ma, James Cheng, and Fan Yu. Dgcl: efficient commu- nication library for distributed gnn training. In Proceedings of the Sixteenth European Conference on Computer Systems, EuroSys '21, pp. 130-144, New York, NY, USA, 2021. Association for Computing Machinery. ISBN 9781450383349, doi: 10. 00.115/3447786.4452333
Michael B Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum. A compositional object-based approach to learning physical dynamics. arXiv, 2016.
Krzysztof Marcin Choromanski. Taming graph kernels with random features. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato. and Jonathan Scarlett (eds.), Proceedings of the 40th International Conference on Machine Learning. volume 202 of Proceedings of Machine Learning Research, pp. 5964-5977. PMLR, July 2023.
Gabriele Corso, Hannes Stark, Stefanie Jegelka, Tommi Jaakkola, and Regina Barzilay, Graph neural networks. Nature Reviews Methods Primers, 4(1):1-13, March 2024. ISSN 2662-8449. doi: 10. 10.10388/43388--24-00222
11

---

Published as a conference paper at ICLR 2026
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2019.
Chinmay Datar, Taniya Kapoor, Abhishek Chandra, Qing Sun, Iryna Burak, Erik Lien Bolager, Anna Veselovska, Massimo Fornasier, and Felix Dietrich. Solving partial differential equations with sampled neural networks. arXiv preprint arXiv:2405.208 2024.
Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and Zico Kolter, End- to-end differentiable physics for learning and control. Advances in neural information processing systems, 31, 2018
Shaan A Desai, Marios Mattheakis, David Sondak, Pavlos Protopapas, and Stephen J Roberts. Port- hamiltonian neural networks for learning explicit time-dependent dynamical systems. Physical Review E, 104(3):034312, 2021.
Gage DeZoort, Peter Battaglia, Catherine Biscarat, and Jean-Roch Vlimant. Graph neural networks at the large hadron collider. Nature Reviews Physics, 5(5):281-303, 2023.
Eva Dierkes, Christian Offen, Sina Ober-Blöbaum, and Kathrin Flaškamp. Hamiltonian neural networks with automatic symmetry detection. Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(6), 2023.
Timothy Dozat. Incorporating Nesterov momentum into Adam. Proceedings of the 4th International Conference on Learning Representations, Workshop Track, May 2016
Weitao Du, He Zhang, Yuanqi Du, Qi Meng, Wei Chen, Nanning Zheng, Bin Shao, and Tie- Yan Liu. Se(3) equivariant graph neural networks with complete local frames. In International Conference on Machine Learning, pp. 5583-5608. PMLR, 2022.
John Duchi, Elad Hazan, and Yoram Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12(61):2111-2159 2011.
Gianluca Fabiani. Random Projection Neural Networks of Best Approximation: Convergence theory and practical applications, February 2024.
Gianluca Fabiani, Francesco Calabrò, Lucia Russo, and Constantinos Siettos. Numerical solution and bifurcation analysis of nonlinear partial differential equations with extreme learning machines. Journal of Scientific Computing, 89(2):44, November 2021. ISSN 0885-7474, 1573-7691. doi:
10.11077/000115-22--1666555
Gianluca Fabiani, Ioannis G. Kevrekidis, Constantinos Siettos, and Athanasios N. Yannacopoulos. RandONets: Shallow networks with random projections for learning linear and nonlinear operators. Journal of Computational Physics, 520: 13433, January 2025. ISSN 00219991. doi: 10. 101 jep.2024.1 1 13433.
David Chin-Lung Fong and Michael Saunders. Lsmr: An iterative algorithm for sparse least-squares problems, SIAM Journal on Scientific Computing, 33(5):2950-2971, 2011.
Evangelos Galaris, Gianluca Fabiani, Ioannis Gallos, Ioannis Kevrekidis, and Constantinos Siettos. Numerical Bifurcation Analysis of PDEs From Lattice Boltzmann Model Simulations: A Parsi- monious Machine Learning Approach. Journal of Scientific Computing, 92(2):34, August 2022. ISSN 7474, 1573-7691. doi: 10. 10.1007/s10915-000 022-01883-y.
Claudio Gallicchio and Alessio Micheli. Graph echo state networks. In The 2010 international joint conference on neural networks (IJCNN), pp. 1-8. IEEE, 2010.
Claudio Gallicchio and Alessio Micheli. Fast and deep graph neural networks. In Proceedings of the AAAI conference on artificial intelligence, volume 34, pp. 1898-3905, 2020.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In Doina Precup and Yee Whye Teh (eds.), Proceedings of the 34th International Conference on Machine Learning, volume 70 Proceedings of Machine Learning Research, pp. 1263-1272. PMLR, August 2017.
12

---

Published as a conference paper at ICLR 2026
Robert Mansel Gower, Nicolas Loizou, Xun Qian, Alibek Sailanbayev, Egor Shulgin, and Peter Richtarik. Sgd: General analysis and improved rates, 2019.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In H. Wallach, H. Larochelle, A. Beygelzimer, d'Alché-Buee E. Fox, and R. Gamett (eds.), Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Anthony Gruber, Kookjin Lee, Haksoo Lim, Noseong Park, and Nathaniel Trask. Efficiently parameterized neural metriplectic systems. In The Thirteenth International Conference Learning Representations, 2025.
Vipul Gupta, Xin Chen, Ruoyun Huang, Fanlong Meng, Jianjun Chen, and Yujun Yan. Graphscale: A framework to enable machine learning over billion-node graphs. In Proceedings of the 33rd ACM International Conference Information and Knowledge Management, CIKM '24, pp. 4514-4521, New York, NY, USA, 2024. Association for Computing Machinery. ISBN 9798400704369. doi: 10.1 00,1144//6677333360021.
Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric numerical integration illustrated by the veörm method. Acta numerica, 12:399 4509, 2003.
William Rowan Hamilton, On a general method in dynamics. Philosophical Transactions of the Royal Society, 124:247-308, 1834.
William Rowan Hamilton, Second essay on a general method in dynamics. Philosophical Transactions of the Royal Society, 144, 1835.
Mohammad Hashemi, Shengbo Gong, Juntong Ni, Wenqi Fan, B. Aditya Prakash, and Wei Jin. A comprehensive survey graph reduction: Sparsification, coarsening, and condensation, 2024.
Quercus Hernández, Alberto Badías, David González, Francisco Chinesta, and Elías Cueto. Structure- preserving neural networks. Journal of Computational Physics, 426: 109950, 2021.
Guang-Bin Huang, Qin-Yu Yu Zhu, and Chee-Kheong Siew. Extreme learning machine: new learning scheme of feedforward neural networks. In 2004 IEEE international joint conference neural networks (IEEE Cat. No. 04CH37541), volume 2. pp. 985-990, 2004.
Guang-Bin Huang, Lei Chen, and Chee Siew. Universal approximation using incremental constructive feedforward networks with random hidden nodes. IEEE transactions on neural networks / publication the IEEE Neural Networks Council, 17:879-92, 2006.
Wei Jin, Xianfeng Tang, Haoming Jiang, Zheng Li, Danqing Zhang, Jiliang Tang, and Bing Yin. Condensing graphs via one-step gradient matching. In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD '22, pp. 720-730, New York, NY, USA, 2022a. Association for Computing Machinery. ISBN 978145039388 3534678.3539429.
doi: 10.1145/
Wei Jin, Lingxiao Zhao, Shichang Zhang, Yozen Liu, Jiliang Tang, and Neil Shah. Graph condensation for graph neural networks. In International Conference on Learning Representations, 2022b. URL
.nt forum htte:/oooommomooox3xxeee
William B. Johnson and Joram Lindenstrauss. Extensions of Lipschitz mappings into a Hilbert space. In Richard Beals, Anatole Beck, Alexandra Bellow, and Arshag Hajian (eds.), Contemporary Mathematics, volume 26, pp. 189-206. American Mathematical Society, Providence, Rhode Island, 1984. ISBN 978-0-8218-50300- 978-0-8218-761 1-4. doi: 10. 1..1990/com/002677774000.
Tim Kaler, Nickolas Stathas, Anne Ouyang, Alexandros-Strin
Iliopoulos, Tao Schardl, Charles Leiserson, and Jie Chen. Accelerating training and inference of graph neural networks with fast sampling and pipelining. In D. Mareulescu, Y. Chi, and C. Wu (eds.), Proceedings of Machine Learning and Systems, volume 4, pp. 172-189, 2022.
D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization, International Conference Learning Representations ICLR 2015, 2015.
13

---

Published as a conference paper at ICLR 2026
Oyku Deniz Kose and Yanning Shen. Fast&fair: Training acceleration and bias mitigation for
GNNs. Transactions on Machine Learning Research, 2023. ISSN 2835-8856. URL https
//openeeveew.eer net/ Eor.mmid- /
Manoj Kumar, Anurag Sharma, Shashwat Saxena, and Sandeep Kumar. Featured graph coarsening with similarity guarantees. In Andreas Krause, Emma Brunskill, Kyunghyun Cho. Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett (eds.), Proceedings of the 40th International Conference Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 17953-17975. PMLR, July 2023.
Kookjin Lee. Nathaniel Trask, and Panos Stinis. Machine learning structure preserving brackets for forecasting irreversible processes. Advances in Neural Information Processing Systems, 34: 5706,5707,,0 2021.
Ho Chun Leung, Chi Sing Leung, and Eric Wing Ming Wong. Fault and noise tolerance in the incremental extreme learning machine. IEEE Access, 7:155171-155183, 2019.
Siye Li, Zhensheng Sun, Yujie Zhu, and Chi Zhang. Physics- constrasconttr and flow-field -ieed--messsgeee informed graph neural network for solving unsteady compressible flows. Physics of Fluids, 36(4), 2024.
Zhiqi Lin, Cheng Li, Youshan Miao, Yunxin Liu, and Yinlong Xu. Pagraph: Scaling gon training on large graphs via computation-aware caching. In Proceedings of the 11th ACM Symposium on Cloud Computing. SoCC 20. pp. 401-415, New York, NY, USA, 2020. Association for Computing Machinery. ISBN 9781450381376 doi: 10.1145/341911111 19111.3321112181
Dong C Liu and Jorge Nocedal. On the limited memory bfigs method for large scale optimization. Mathematical programming, 45(1):503-528, 1989.
Liyuan Liu, Haoming Jiang, Pengcheng He, Weizhu Chen, Xiaodong Liu, Jianfeng Gao, and Jiawei Han. On the variance of the adaptive learning rate and beyond, 2021.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization, 2019.
Michael Lutter, Christian Ritter, and Jan Peters. Deep lagrangian networks: Using physics as model prior for deep learning. In International Conference on Learning Representations, 2019.
Antonio Marino, Claudio Pacchierotti, and Paolo Robuffo Giordano. A gated graph neural network approach to fast-convergent dynamic average estimation. ACM Trans. Intell. Syst. Technol,, March 2025. ISSN 2157-6904. doi: 10. 145/322585. Just Accepted.
Xiangrui Meng, Michael A Saunders, and Michael W Mahoney. Lsmn: A parallel iterative solver for strongly over-or underdetermined systems. SIAM Journal Scientific Computing, 36(2): C95-C118, 2014a.
Xiangrui Meng, Michael A. Saunders, and Michael W. Mahoney. LSRN: A parallel iterative solver for strongly over- or underdetermined systems. SIAM Journal on Scientific Computing, 36(2): C95-C118 18, January 2014b. ISSN 1064-8275, 1095-7197. doi: 10. 11117722200668800
Philip M Morse. Diatomic molecules according the wave mechanies. ii vibrational levels. Physical review, 34(1):57, 1929.
Mohammad Amin Nabian, Chang Liu, Rishikesh Ranade, and Sanjay Choudhry. X-meshgraphnet: Scalable multi-scale graph neural networks for physics simulation, 2024.
Amrit Nagarajan and Anand Raghunathan. FASTRAIN-GNN: Fast and accurate self-training for graph neural networks. Transactions on Machine Learning Research, 2023. ISSN 2835-8856.
Sina -bor-Bloe and Christian Offen. Variational learning of Euler Lagrange dynamics from data. Journal of Computational and Applied Mathematics, 421: 14780. 2023.
C. Offen and S. Interdisciplinary Journal of Nonlinear Science, 32(1):013122, 2022.
Symplectic integration of learned Hamiltonian systems. Chaos: An
14

---

Published as a conference paper at ICLR 2026
Christopher C Paige and Michael A Saunders. Lsqr: An algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software (TOMS), 8(1):43-71, 1982.
Y-H Pao and Yoshiyasu Takefuji. Functional- link net computing: theory, system architecture, and functionalities. Computer, 25(5):76-3) 79, 1992.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary Vetito Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems 32, pp. Associates, Inc., 2019.
Curran
Jiang-Zhou Peng, Yue Hua, Yu-Bai Li, Zhi-Hua Chen, Wei-Tao Wu, and Nadine Aubry. Physics- informed graph convolutional neural network for modeling fluid flow and heat convection. Physics of Fluids, 35(8), 2023.
Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter Battaglia. Learning mesh-based simulation with graph networks. In International Conference on Learning Representations, 2021.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007.
Ali Rahimi and Benjamin Recht. Uniform approximation of functions with random bases. In 2008 46th annual allerton conference on communication, control, and computing, pp. 555-561. IEEE, 2008.
Atamert Rahma, Chinmay Datar, and Felix Dietrich. Training Hamiltonian neural networks without backpropagation. In Workshop on Machine Learning and the Physical Sciences. NeurlPS 2024, November 2024.
Isaac Reid, Krzysztof Choromanski, Eli Berger, and Adrian Weller. General graph random features, 2023a.
Isaac Reid, Krzysztof M Choromanski, and Adrian Weller. Quasi-monte carlo graph random features. In A. Oh, Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Neural Information Processing Systems, volume 36, pp. 17776-1779, Curran Associates, Inc., 2023b.
Martin Riedmiller and Heinrich Braun. A direct adaptive method for faster backpropagation learning: The rprop algorithm. In IEEE international conference on neural networks, pp. 586-591. IEEE, 1993.
Herbert E. Robbins. A stochastic approximation method. Annals of Mathematical Statistics, 22: 400-407, 1951.
Frank Rosenblatt. Perceptions and the theory of brain mechanisms. Spartan books, 1962.
Fabian J. Roth, Dominik K. Klein, Maximilian Kannapinn, Jan Peters, and Oliver Weeger. Stable Port-Hamiltonian Neural Networks, February 2025.
Alvaro Sanchez-Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian Graph Networks with ODE Integrators, September 2019.
Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovee, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International conference PMLR, 2020.
machine learning. pp.
Michael Schlichtkrull, Thomas N. Kipf, Peter Bloem, Rianne van den Berg, Ivan Titov, and Max Welling. Modeling relational data with graph convolutional networks. In The Semantic Web: 15th International Conference, ESWC 2018, Heraklion, Crete, Greece, June 3-7, 2018, Proceedings, pp. 593-607, Berlin, Heidelberg, 2018. Keringer-Verlag ISBN 978-3-319.93466-
15

---

Published as a conference paper at ICLR 2026
Robin M Schmidt, Frank Schneider, and Philipp Hennig. Descending through crowded valley benchmarking deep learning optimizers. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 9367-9376. PMLR, 2021.
Wouter F Schmidt, Martin A Kraaijveld, Robert PW Duin, et al. Feed forward neural networks with random weights. International conference on pattern recognition, pp. 1-1. IEEE Computer Society Press, 1992.
Yingxia Shao, Hongzheng Li, Xizhi Gu, Hongbo Yin, Yawen Li, Xupeng Miao, Wentao Zhang, Bin Cui, and Lei Chen. Distributed graph neural network training: A survey. ACM Comput. Surv., 56 (8), April 2024. ISSN 0360-0300. doi: 10. 10.1145/3648358
Vinay Sharma and Olga Fink. Dynami-cal graphnet: A physics-informed graph neural network conserving linear and angular momentum for dynamical systems, 2025.
Noam Shazeer and Mitchell Stern. Adafactor: Adaptive learning rates with sublinear memory cost. PMLR, 2018.
In International Conference on Machine Learning, pp.
Khemraj Shukla, Mengjia Xu, Nathaniel Trask, and George E. Karniadakis. Scalable algorithms for physics-informed neural and graph networks. Data-Centric Engineering, 3:e24, 2022. doi: 10.1017/dce.2022.24
lya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. On the importance of initialization and momentum in deep learning. In International conference on machine learning, pp. 1139-1147. PMLR, 2013.
Joshua Tenenbaum, Vin de Silva, and John Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500) 2000.
Abishek Thangamuthu, Gunjan Kumar, Suresh Bishnoi, Ravinder Bhattoo, NM Krishnan, and Sayan Ranu. Unravelling the performance of physics-informed graph neural networks for dynamical systems. Advances in Neural Information Processing Systems, 35:3691-3702, 2022.
Tijmen Tieleman and Geoffrey Hinton. Rmsprop: Divide the gradient by a running average of its recent magnitude. lecture 6.5, 2012.
Alicia Tierz, Icfar Alfaro, David González, Francisco Chinesta, and Elías Cueto. Graph neural networks informed locally by thermodynamics. Engineering Applications of Artificial Intelligence, 144:110108 March 2025. ISSN 09521976. doi: e...006/..enaapai...a 2025. 10108.
Alan John Varghese, Zhen Zhang, and George Em Karniadakis. SympGNNs: Symplectic Graph Neural Networks for identifying high-dimensional Hamiltonian systems and node classification. Neural Networks, 187:107397, July 2025. ISSN 08936080. doi: 101106/j.nuunet2200 107397.
Clément Vignac, Andreas Loukas, and Pascal Frossard. Building powerful and equivariant graph neural networks with structural message-passing. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 14143-14155. Curran Associates, Inc., 2020.
Xinchen Wan, Kaiqiang Xu, Xudong Liao, Yilun Jin, Kai Chen, and Xin Jin. Scalable and efficient full-graph gnn training for large graphs. Proc. ACM Manag, Data, 1(2), June 2023. doi: 10. 1145/ 3589288.
Lei Wang, Qiang Yin, Chao Tian, Jianbang Yang, Rong Chen, Wenyuan Yu, Zihang Yao, and Jingren Zhou. Flexgraph: flexible and efficient distributed framework for gnn training. In Proceedings of the Sixteenth European Conference on Computer Systems, EuroSys '21, pp. 67-82, New York, NY, USA, 2021. Association for Computing Machinery. ISBN 978445038349,, doi: 10.1 10.1145/3447786.345556
Shaocong Wang. Yi Li, Dingchen Wang, Woyu Zhang, Xi Chen, Danian Dong, Songqi Wang, Xumeng Zhang, Peng Lin, Claudio Gallicchio, et al. Echo state graph neural networks with analogue random resistive memory arrays. Nature Machine Intelligence, 5(2):104-113, 2023.
16

---

Published as a conference paper at ICLR 2026
Nicholas Watters, Daniel Zoran, Theophane Weber, Peter Battaglia, Razvan Pascanu, and Andrea Tacchetti. Visual interaction networks: Learning a physics simulator from video. Advances in neural information processing systems, 30, 2017.
Shiying Xiong, Yunjin Tong, Xingzhe He, Shuqi Yang. Cheng Yang. and Bo Zhu. Nonseparable symplectic neural networks. In International Conference on Learning Representations, 2021.
Tianju Xue, Zhengtao Gan, Shuheng Liao, and Jian Cao. Physics-embedded graph network for accelerating phase-field field simulation of microstructure evolution in additive manufacturing. npj Computational Materials, 8(1):201, 2022.
Jianbang Yang, Dahai Tang, Xiaoniu Song, Lei Wang, Qiang Yin, Rong Chen, Wenyuan Yu, and Jingren Zhou. Gnnlab: a factored system for sample-based gnn training over gpus. In Proceedings of the Seventeenth European Conference on Computer Systems, EuroSys '22, pp. 417-434, New York, NY, USA, 2022. Association for Computing Machinery. ISBN 978 9781450391627, doi: 10.1 10.1145//99322..3119555
Youn- Yool Yu, Jeongwhan Choi, Jaehyeon Park, Kookjin Lee, and Noseong Park. Piorf: Physics- informed ollivier-ricci flow for long-range interactions in mesh graph neural networks, 2025.
Matthew D Zeiler. Adadelta: adaptive learning rate method. arXiv preprint arXiv:1212.5701. 2012.
Lizhi Zhang, Zhiquan Lai, Shengwei Li, Yu Tang, Feng Liu, and Dongsheng Li. 2pgraph: Acceler- ating gnn training large graphs gpu clusters. In 202/ IEEE International Conference on Cluster Computing (CLUSTER), pp. 103-113, 2021. doi: 10.1109/Cluseer44895522099 1.00036.
Rui Zhang, Yuan Lan. Guang-Big -aing--Binne Huang, and Zong-Ben Xu. Universal approximation of extreme learning machine with adaptive growth of hidden nodes. IEEE Transactions on Neural Networks and Learning Systems, 23(2):365-371, 2012.
Xin Zhang, Yanyan Shen, and Lei Chen. Feature-oriented sampling for fast and scalable gnn training. In 2022 IEEE International Conference on Data Mining (ICDM), pp. 723-732, 2022a. doi: 10.1 100/CCDM54844.22222200
Zhen Zhang, Yeonjong Shin, and George Em Karniadakis. GFINNs: GENERIC formalism informed
neural networks for deterministic and stochastic dynamical systems. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 38002229):20210007, August 2022b. ISSN 1364-503X, 1471-2962. doi: 10. 20998/staa2021.0000 1.0077.
Xinge Zhao and Chien Chern Cheah. Ensuring reliable learning in graph convolutional networks: Convergence analysis and training methodology. IEEE Transactions on Artificial Intelligence, pp-
1-15, 2025. doi: 10.1 10.110/AA.2225555555555
Yingxue Zhao, Haoran Li, Haosu Zhou, Hamid Reza Attar, Tobias Pfaff, and Nan Li. A review of graph neural network applications in mechanics-related domains. Artificial Intelligence Review, 57 (11):315, 2024.
Hongkuan Zhou, Da Zheng, Israt Nisa, Vasileios Ioannidis, Xiang Song, and George Karypis. Tgl: general framework for temporal gnn training billion-scale graphs. Proc. VLDB Endow, 15(8): 1572-1580, April 2022. ISSN 2150-8097,
doi: 10.14778152933.352994..
Zeyu Zhu, Peisong Wang, Qinghao Hu, Gang Li, Xiaoyao Liang, and Jian Cheng. Fastgl: A gpu- efficient framework for accelerating sampling- based gnn training large scale. In Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume ASPLOS '24, pp. 94-110, New York, NY, USA, 2025. Association for Computing Machinery. ISBN 979840070391 1. doi: 10.1 1e.1111//222881366777
Nikita Zozoulenko, Thomas Cass, and Lukas Gonon. Random feature representation boosting, In Forty-second International Conference on Machine Learning, 2025.
17

---

Published as a conference paper at ICLR 2026
Table A.5: Notation used in Section 3 and in the appendix.
Problem setup
N d M
Number of nodes/particles in the system
Spatial dimension
Number of training data points
Ne
Number of edges in the graph
4.P E TRd-N
Positions and momenta of all the nodes
Pi E Rd
Positions and momenta of the eth node True Hamiltonian function
R2d-N
H: R2d-N
Predicted Hamiltonian function Symmetric adjacency matrix
A
RF-HGN setup
V E
Set of node feature encodings Set of edge feature encodings Number of node features
dy E N
dr E N
Number of edge features Latent (hidden) dimension for encoding node and edge features Latent (hidden) dimension for encoding messages Size of the input to the linear layer (here also the network width) Invariant representation Node features for the hth node
dn N
dM E
dp EN
h. E Rd
1/2 E R
6gj E Rd
Edge features for the edge (i,j) with i Set of source nodes connected to destination node dense layer for encoding the node features dense layer for encoding the edge features Message encoder Weights and biases of the node encoder Weights and biases of the edge encoder
N
Øv Rd + Rd.
OF Rd 7 Rd
0M 2d, + Rd
Wy Rdaxds by E Reda
WE Rdaxds by E Rd
WM E Rdd XR bM E TRdar
Weights and biases of the message encoder Weights and biases of the linear layer aggregated incoming message of node Node feature encoding for the hh node
WL Rdd by E
m E Rdu
h E Rd
EE E IRdh
Edge feature encoding for the edge defined by nodes i,
MM E
Rdu
Constructed message encoding
bG E Rdc
Global encoding of the graph
APPENDIX
A NOTATION
We define the notation we use in the main text and in the Appendix in Table A.5
B ADDITIONAL DETAILS ON RANDOM-FEATURE HAMILTONIAN GRAPH NETWORKS
We discuss the problem setup used in Figure i Appendix B.1 how to construct rotation-invarian representation for spatial dimension d = Appendix B.2, the algorithm for forward pass for RF-HGN in Appendix B.3 and run-time and memory complexity in Appendix B.4
B. ADDITIONAL DETAILS FOR FIGURE
The initial conditions for training the networks in Figure 2 are generated by displacing the positions by dg U--0.5, +0.5) and with momenta p U(-2, +2) from some fixed reference frame illustrated
18

---

Published as a conference paper at ICLR 2026
on the left of Figure 2 The system is integrated using symplectic Störmer- Verlet Hairer et al. (2003)
with At = 10 for 100 steps using the true dynamics and predictions of invariant and non-invarian variations of (SWIM) RF-HGN.
B.2 ROTATION-INVARIAL REPRESENTATION FOR SPATIAL DIMENSION d ===
The details on obtaining iattionninvariant representations for the spatial dimension d are discussed in Section 3
Here, we extend this approach to the spatial dimension d = 3, leveraging the classical Gram-Schman orthogonalization method.
We pick E R and 92 E R$ the reference vectors and let
1i and eg Teall is the first basis vector of the new frame. If lefég] > (near colinear), set ely x el (cross product). We then W22 (Gram-Schmidt) to compute the second basis vector wall
project u2 E/ proje, 5/ and scale e2
62, where proje, b - (el E1)ei is the projection of e/2 onto E₁. We then construct E3 = E1 X 62, where the cross-product. Finally, we define the orthonormal basis B [e] 62 Ez]. We set 0.98 if not specified otherwise.
One can uniquely identify the first two points, independent of node ordering or orientation, as the ones closest to the mean g. In case of ties, we select the point with the smallest angle relative to the first coordinate axis centered at q. If the ties persist, we can then select the point with the smallest angle relative to the second coordinate axis.
B.3 ALGORITHM
Here, we outline the algorithm for the forward pass of the Random-Feature Hamiltonian Graph Network (RF-HGN) using the notation introduced in Section 3
The parameters of all dense layers Φv, &gg, dı are computed leveraging random sampling techniques and last layer parameters WL and by are computed using least squares (see Section (.3). We denote the set of neighbors that transmit information to node j by Ng. In the following, we use a single subseript, for instance, for U; to double subscript, for instance, for eig, to denote that we compute eiy for i, E (1,..., N} and > j. and set Eji Cij-
Algorithm B.1: Forward pass for RF-HGN:
denote that compute v; for all values of i 21,2.... N) for brevity. Also,
Input: Positions and momenta of the bodies in spatial dimension (p, q E 1R2 adjacency matrix A E x
Output: Approximation
Hamiltonian A
Parameters: Node/edge encoder dimension da E N and message encoder dimension das N
Ts. p. E 2d +- encode wwwccodeiinvrrianeeeee q) for each E {l...., .N} (Encode tramslation- and rotation-invariance)
0, - [4, ]:]
E 2
(Node features)
Cgj [(ā. 4 |g,
Till] E Rd+1
(Edge features)
h (v: E Rd
[Node encoding
EE de(en) E Rdh
(Edge encoding
MM + Φm (I(N))T (1E)71)) Iddar
h
[Message encoding]
my + EieN, hi e Rda
[Message passing (local pooling)
hG Σ3= [(h)) (m,)TTT Rdi, where dr. dr + dM (Message passing (global
pooling))
H. + Wr. ha + bi, where W₁, Idd and
IR
[Linear layer)
return H
The forward pass discussed here is independent of how the network parameters are computed. The training leverages random sampling, automatic differentiation to compute gradients of the Hamiltonian with respect to inputs to compute VH using PyTorch (Paszke al., 2019), and least
19

---

Published as a conference paper at ICLR 2026
squares solvers as described in Section
3.3 The network can be formulated more compactly as
R(G)) - WL Σ
Σ VEN
ΦA (0<(o((ss) + PL.
(B.7)
BYEV
where || is concatenation, given graph G = (V, with node feature set V and edge feature set E as explained in Section problem setup.
B.4 RUN-TIME AND MEMORY COMPLEXITY OF TRAINING
We use the notation defined in Table
The bottleneck of the run-time complexity is described in Section 3.3 Encoding the translational symmetry requires mean-shift of the particles which can be done in ((MNd) because for each system first the mean value of the positions has be computed and then the values are updated which can all be done linearly in M, N, and d given d positions we have to shift. For encoding rotational symmetry we have implemented Gram-Schmin orthogonalization, which is in (MMNd))
Run-time complexity:
Also note that there are M(M-1) pairs of data points choose from when sampling random features with SWIM. In practice, we do not consider all possible pairs, but rather subsample this set uniformly
by choosing the candidate number of pairs to much less than the theoretically possible number of pairs, and still results in robust sampling method.
[W]
M, where [W] is the number of neurons. This
M
Memory complexity: Memory requirements for training set of size M graphs include O(d-N (id-N-JJM -M) node features, O(d-N; . .M) for edge features, and ( for the sparse adjacency matrix, assuming For sparsity, we assume O(1) number
the graph stays the same for each example in the training of neighbors for each node. The three dense layers (node, edge, and message encoders incur costs of (dn-dv)) Odd -(), O(d 1ala--dp, dp.). The linear readout layer adds further ((z)) Odh + +dM) (dd))
Unlike gradient-descent-basee iterative optimization schemes, we only need compute the gradients of the Hamiltonian H with respect to inputs, and not with respect to parameters. For this, we additionally need to store the partial derivatives of the output with respect to the input of each dense layer for bocckpropagation. This amounts to an additional cost of O(dy d M) for the partial derivatives of the global graph value with respect to inputs.
For fixed spatial dimension dataset size M, the number of nodes N. and the number of edges in graph Ne. the total memory footprint during training is O(M + Ne)).
and network width d,, since the dominant terms depend on
). If we further assume zero-shot generalization with fixed training system, then the total memory requirement is in O(M) and the geometry of the system (the number of nodes N and edges N.) can grow independently of this training.
C DATASETS
Table C.6lliss lists summary information of the datasets used in our experiments, which are explained in more detail in the following subsections. All the constants (masses, spring, and Lennard-Jone constants) are set to one all the experiments, and for the chain and lattice examples we have used relative distances (all positions given displacements relative to the equilibrated state). More information can found the code repository: https: /qiilab com/ fdd Edwresearc//h swi mhgn
20

---

Published as a conference paper at ICLR 2026
Table c.6: Summary of the datasets used in our main experiments. Low and high specify the uniform distribution used to sample the dataset. In Appendix C.1 Appendix and Appendix C.3 we give more details on how we generate the datasets.
Experiment
Train points
Test points
Low
High
Table Tablo
1000 N 6
1000 N G
-0.5
+0..5
3000 N
3000 4
-1.0
+1.0
Figure
NN xy train)
1000 N₂ N,
1000 N₂ N2.
-0.5
+0.55
Figure
yy y test)
1000 N₂ Ny
-0.5
+0.5
Figure
0.16
and
train N)
2000 N 44
2000 4
-1.0
+1.0
Figure
F1.16
2000 N 4
-1.0
+1.0
Table
Figure H.18
810
N
+0.1
and H. 19
Higure
11.20
H.21
and
H121
540
60 N
0.11
+0.1
Section 4.4
10000 N -
100 N
C.I BENCHMARKING AGAINST SOTA OPTIMIZERS
The target Hamiltonian is
Nu
Ny
1
(g =
2
(C.8)
Nz
Ny-1
Ny Nz-1
B 16g41
2i5ll +
Σ
Σ
B 1 -
where Pty,Pi R3, and aij. Pú E R denote masses and spring constants, respectively. All constants equal to one if not specified otherwise. N₂ and Ny set to three to build 3x3 lattice structure (with number of total nodes N = 9), which moves in 3D (d = 3). We generate a synthetic dataset of 2000 structures (graphs) with their true time derivatives (qi, Ps: dt. P: i=1 where q. Pi. 4i P. 1dd-N Vi. We first set all gi P: to be in the equilibrium state. Then displacements dq, and dp; from the uniform distribution U(-0.5, +0.5), and compute q < qi + dq; and p; <- P: + dp,. We then compute the ground truths q: Ps using Equation (1) and the ground truth gradient VHi We shuffle and split the dataset into train (1000) and test (1000) sets. All the errors reported in Table are the average test errors of three independent runs using different seeds. The total number of training and test points then becomes 1000 N W.dd. d 54000 each.
sample the
Additional to the standard spring potential V(r) 1rr2² given distance with spring constant / we use an anrharmonic spring potential V(r) = _,,r² + 1m with nonlinearity coefficient r) = and the Morse potential (r) D(1 exp(-ar))2 (Morse 1929) with well-depth 1 and potential-width
1.0 for the 2D chain potential experiments in Table 3 (also in Figure 10 H.11
H.12 H.13 H.14 and 15). The data generation follows the same procedure explained above, with N 5, resulting in a total of 3000 N d 2 60000 train and test points each.
C.2 ZERO-SHOT GENERALIZATION AND COMPARISON RANDOM FEATURE METHODS
The experiment in Figure 5 uses the same procedure explained in Appendix C.I with N₂ and N, set to two, three, and four to build 2x2, 3x3, and 4x4 lattice structures.
For the experiment in Figure H6 the procedure is again similar, but the structure of the experiment and data is different (an open chain). The target function for the open chain system given in Section 4.2 as
Pp 12
Hzlg.p)
Σ
B.11 9.ll7).
(C.9)
where qs Pi E RR, a; B5 E R are positions, momenta, masses, and spring constants in the system, respectively, for i E NJ. All constants equal one if not specified otherwise. H2. N scaled from exponentially from 2 to 22 in the experiment, which always moves in 2D (d = 2). For each N, we again generate synthetic dataset of 4000 structures (graphs) with their true time
21

---

Published as a conference paper at ICLR 2026
Table D.7: Model parameters (see Figure
used in Section 4.1
Model
Encoder width (d,,)
Network width (d,))
Activation
Precision
HGN (fg
48
384
softplus
single
Table D.8: Hyperparameters used in Section 4.1 and Section 4.2 are listed for SWIM. Dri ver and reond ((2 reg.) are the parameters of t torch... linallg istag PPaszke et al. 2019). Resample
duplicates specifies to resample till we get a unique pair of points in the SWIM algorithm (Bolager
ct al. 2023
Optimizer
Driver
Parameter sampler
Resample duplicates
reg.
SWIM SWlMg al. 2023
GELS
rolu
Trie
le.6
derivatives {gi: Pi, fb }t 2000 where qi,Pi, qs, P. Vi. We first set all Ge, to be in the equilibrium state. Then sample the displacements dq; and dp, from the uniform distribution U--1.0, +1.0). and compute & + dq; and p; + + dp;. We then compute the ground truths q; p; using Equation (I) and the ground truth gradient vH2. We shuffle and split the dataset into train (2000) and test (2000) sets.
For the molecular dynamics scenarios with the Lennard-Jones (LJ) potential the Hamiltonian is defined as
Ipsll
Hs(g,p)
V (11 %.l)),
(C.10)
ag i=1
j=i+1
where
12
Vi
4e
and where raj = 11 1:. We set the parameters ag 6, a 1.0 and the cutoff to 2.0 when computing the dynamic edge indices (Figure H. H.21 and static edge indices when training with 9 particles and testing with particles (Figure H.18| H. Table H.32 H.33). For the static edge experiment, we generated 300 trajectories with 9 particles, and for the dynamic edge experiment, we generated 200 trajectories with 36 particles with the & displacement specified in Table|C.6 state with momenta set to zero. Each trajectory is simulated for 50 time steps with A/ = 5 10-3 and snapshots are taken every 20th step with a train-test ratio of 0.9.
the equilibrium
C.3 BENCHMARKING AGAINST SOTA ARCHITECTURES
To benchmark model, we considered the N-body spring system from Thangamuthu al. (2022) for which details are available the original work. Nonetheless, we mention the key properties of the dataset for completeness.
A system of N bodies with equal masses, connected by elastic springs such that each body has two connections and the system forms closed loop. The system's 's physical behavior additionally depends on the spring's 's stiffness and its undeformed length, both are set to one. Initial positions go are sampled as qo (0, 2) and initial momenta Po are sampled as po (O, 0.1) and subsequently mean-centered. The symplectic Störmer- (Hairer et al. 2003) integrator with a timestep of 10-3 is used to generate 100000 datapoints, which are subsampled to 100 datapoints. The approach is repeated for 100 trajectories obtain dataset that is split in 75:25 ratio for training and validation set. Unlike the original work, the test data we use consists of only one trajectory because with 100 trajectories, we were often experiencing failed simulations with the existing Adam-trained benchmarks (in particular with the LGN architecture), which significantly hinders comparison with method.
D TRAINING AND ZERO-SHOT INTEGRATION SETUP
D.1 BENCHMARKING AGAINST SOTA OPTIMIZERS
Table D.7 and Table D& list the model and SWIM hyperparameters, respectively. Table (.9 the hyperparameters used for the SOTA optimizers listed in Table All the optimizers run once
22

---

Published as a conference paper at ICLR 2026
Table D.9: Hyperparameters used in Section 4.1 and Section 4.2 are listed for SOTA optimizers. SGD(+m) represents both SGD Robbins 1951 I and SGD+momentum (Sutskever et al. 2013). The momentum parameter is set to 0.9. Avg. SGD specified the Averaged SGD (Gower al. 2019). Default values are given in the first row (Defaults) for all the optimizers not present in this table, but are listed in Table 1 #steps is the number of total iterations (one iteration per batch). If LR schedule is specified, exponential decay is used as the learning rate scheduler. All optimizers use the kaiming hasmiinggoorralllm normal (Paszke et al., 2019) weight initialization. 1 regularization (?² reg.) is specified
using the weight _eccttd_____y__ee parameter (Paszke al. 2019). Full batch size is 1000.
Optimizer
#steps
Batch size
LR schedule
LR (sturt, end)
1 reg.
Defaults
10000
256
Yes
le 5e
le-6
SGD(+m))
10000
256
Yes
5e-4, 5e-5
le-6
Avg. SGD
10000
256
Yos
5c-4, 5c-5
0
Adadelta
10000
256
No
le-1 (lixed)
le-6
Rprop
2560
Full
No
1e-2 (fixed)
0
RMSprop 10000
256
Yes
lc-2, Se-5
0
LBFGS
100
Full
No
(fixed) D
Table D.10: Parameters used in Table (also in Figure H.10 is the total number of time steps, At is the time step size.
H.11 H.12
H.13
H.14
H.15
#steps
Parameter
Value
Symplectic solver Störmer-Verlet (Hairer
2003
#steps
Iea
with the default settings that are optimized initially for the Adam optimizer (Kingma & Ba 2015).
and tuned further with multiple iterations. Note that we only want give "time to solution" with similar accuracies in order to compare the SOTA optimizers against our method, since the iterative routines can be done arbitrarily long and may be tuned further reach lower approximation errors than our method with excessive hyperparameter tuning and larger number of iterations for each optimizer-at at the cost of even longer training times. The SGD family (SGD (Robbins) I95I). SGD+momentum (Sutskever et al. 2013). and Averaged SGD (Gower et al. 2019)) required lower learning rate starts than the ceaptive--radieen based optimizers, otherwise they led to NaN (Not a Number) results. Even with a very low learning rate, starting at 5e-4, they all produced one NaN value out of three experiments, which shows their instability and difficulty in setup. In our results, we therefore only average the two valid results of the SGD-family. In Averaged SGD (Gower et al.) 2019), the averaging may have acted as an implicit regularizer, and required no weight decay to perform similarly. Also, regularization LBFGS (Liu & Nocedal 1989), and RMSprop (Tieleman & Hinton 2012). Adadelta (Zeiler 2012) an adaptive method that dynamically scales updates; therefore, it does require any scheduler.
not necessary for Rprop (Riedmiller & Braun 1993).
Also, Rprop (Riedmiller & Braun] 1993) uses the sign of the gradients and adapts the step size dynamically, which makes suitable be used with large batch updates and no scheduler. Since uses full batch updates (with batch size of 1000), its number of gradient steps is reduced to provide around the same epoch as the other optimizers. LBFGS is second-orde method and outperformed the other optimizers with only 100 steps using full batch updates and no learning rate scheduler.
Table D.II D.14 D.12, and D.14 list model hyperparameters for the chain potential experiment in Table 3 for training. Table To the integration hyperparameters used in the same experiment when zero-shot testing. Note that during testing, we apply constant gravitational force [0, 0.075] to every node in the negative y-axis direction.
D.2 ZERO-SHOT GENERAI AND COMPARISON OF RANDOM FRATURE METHODS
Table D.15. Table D.8 and Table 16 experiments in Section
( list the model, SWIM, and ELM hyperparameters used for the respectively. For the zero-shot evaluation presented Figure 3 we have
Table D.11: Model parameters (see Figure 3) used in Table
Model
Encoder width (d.)) Network width (dz) Activation
Precision
HGN t0g.
64
1024
scftplus
deuble
23

---

Published as a conference paper at ICLR 2026
Table D.12: Hyperparameters used in Table are listed for SWIM. Dri vee and rcond (?2 reg.) are the parameters of torch. linalg...... lssslsssg (Paszke et al. 2019). Resample duplicates specifies to resample till we get a unique pair points the SWIM algorithm Bolaagrr et al. 2023).
Optimizer
Driver
Parameter sampler
Resample duplicates
1 reg.
Bolager
2023
GELS
Frue
le-10
Table D.13: Hyperparameters used in Table 3 are listed for ELM (Huang et al., 2004 Driver and rcond (l2 reg.) are the parameters of torch... linalg. . listaq (Paszke al. 2019). Bias low and high specify the uniform distribution of low and high values, from which the biases of the random feature layers sampled. The weights sampled using the standard normal distribution as explained Section 3.3
Optimizer
Driver
Bias low
Bias high
t"
(Huang
GELS
trained (SWIM) RF-HGN ten times with different random seeds (also see Figure H.8.
Figure H.9 and
Table H.3I), and used the pretrained (SWIM) RF-HGN model with the median test error to evaluate on the zero-shot test cases in order to avoid any statistical bias, this is a random feature method.
Table D.17 lists the parameters used to integrate the system in Figure 6 and Figure H.17
For the molecular systems, Table D.19 lists the model parameters, Table D.20 lists SWIM hyperpa- rameters, D.21 ELM hy pyprpaaraetter, and Table|D.22) lists Adam hyperparameters. No early stopping triggered in these experiments. Table D.18 lists the parameters used to integrate all the molecular systems presented in this paper.
D.3 BENCHMARKING AGAINST SOTA ARCHITECTURES
Table D.23 lists the model parameters, Table D.24 lists SWIM hyperparameters, Table D.25 lists
SOTA architecture hyperparameters used in Section 4.4
E
HARDWARE
Table E.26 lists hardware used for all the experiments presented in Section The experiments presented in section 4.1 Figure (trannin of the 2x2 and 3x3 lattice systems), and Section 4.4 are conducted on a CUDA GPU. Figure 5 (training of the 4x4 lattice), Figure (testing,, and Figure 16 conducted CPUs because of their larger memory requirements.
F ABLATION STUDIES
We vary the widths of the encoders and linear layers to understand how they affect the mean squared error defined on the true and predicted solutions. We experimented with the Hamiltonian H2 explained in Section 4..] chain of 8 nodes in 2D. We evaluated the model on the training and testing 2000 samples each in the phase space and report the test errors. The message encoder's width is chosen by subtracting the width of the hidden layer from the width of the linear layer in the ablation study, i.c., dı d, d,, Figure F.7 and Table F.27 reveal that increasing either the linear layer width d, or the hidden dimension d, while keeping the other parameter fixed consistently reduces the mean squared error. In the same experiment, we also computed the condition number k(Z) a, in terms of the singular values > ... > on of the matrix Z associated with the linear system in Equation (6) assess the sensitivity of the solution. We avoided the bias term when computing the condition, as it was only used to fit the integration constant in practice. Table F.88 reveals larger values
with
Table D.14: Hyperparameters used in Table are listed for Adam.
Optimizer
#steps
Batch size
LR schedule
LR (start, end)
reg.
Initialization
Patience
Adam
10000
256
exponential decay
10-2,
1e-6
Kaiming nermal
1000
24

---

Published as a conference paper at ICLR 2026
Table D.15: Model parameters used in Section 4.2 to the last linear layer in both RF-HNN (Bertalan et al. 2019 Greydanus et al. 2019) and RF-HGN.
The network width specifies the size of the imput
Model
Encoder width (da.) Network width
Activation
Precision
RE-HIG (fig.
saftplus soštplus
singla
RE-HNN
512
single
Table 16 Hyperparameters used Section 4.2 are listed for ELM (Huang et al. 2004). Driver
and rcond (1² reg.) are the parameters of torch. linalg. lsttag (Paszke al. 2019). Bias low and high specify the uniform distribution of low and high values, from which the biases of the random feature layers sampled. The weights are sampled using the standard normal distribution
as explained Section 3.3
Optimizer
Driver
Bias low
Bias high
1
Huang etal.
2004
GELS
+1.0
Io-6
Table D.17: Parameters used in Figure 6 (and consequently Figure of time steps, At is the time step size.
H.17
#steps is the total number
Parameter
Value
Symplectic solver Störmer-Verlet (Hairer et al
2003
#steps
5000
At
le-3
Table D.18: Parameters used in Figure H.18 and Figure H. (and consequently Table 2 H.32
H.33
Figure H.19 H.21). #steps is the total number of time steps, At the time step size.
Parameter
Value
S Symplectic solver
Störmer-
Hairer et al
2003
#sleps
At
le-5
Table D.19:
Model parameters used Table 2
H.32
H.33
Figure H.18 H.19 (top row) and in
Figure H.20
H.2 (bottom row) are listed for the RF
RF-HGNN HGNs and Adam-HGN.
Model
Encoder width (da)
Network width
Activation
Precision
HGN
40
800
softplus gelu
single single
HIGN
32
256
Table D.20: Hyperparameters used in Table 2
H.32
H.33
Figure H.18 H.19 and Fig-
H.20 H.21 are listed for SWIM. Driver and reond (12 reg.) the parameters of
terch. inalg. latsq (Paszke et al. 2019). Resample duplicates specifies to resample till get unique pair of points in the SWIM algorithm (Bolager al. 2023)
Optimizer
Driver
Parameter sampler
Resample duplicates
reg.
SWlafBalagerrea et al. 2023
GRT.89
relu
True
le-10
Table D.21: Hyperparameters used Table
H.32
H.33
Figure H.18 H. and in Figure H.20
H.21 are listed for ELM (Huang et al. 2004). Driver and rcond (1a reg.) are the parameters of torch. wwi.nlnni..... ls (Paszke et al. 2019). Bias low and high specify the uniform distribution of low and high values, from which the blases of the random feature layers sampled. The weights are sampled using the standard normal distribution as explained Section
3.3
Optimizer
Driver
Bias low
Bius high
Husng
2004
Table D.22: Hyperparameters used in Table H.2T are listed for Adam.
H.32
H.33
Figure
H.18 H.18
H.19 and in Figure H.20
Optimizer
#steps
Batch size
LR schedule
LR (start, end)
reg.
Initialization
Patience
Adam
10000
exponential decay
Ie-3, Se-5
le-10
Kaiming nommal
500
Table D.23: Random feature Hamiltonian graph network (RF-HGN) parameters used in Section 4.4
The network width specifies the size of the input to the last linear layer in the RF-HGN.
Model
Encoder width (dx)
Network width
Activation
Precision
RI-HO HGN (lg.
32
512
soštplus
double
25

---

Published as a conference paper at ICLR 2026
Table D.24: Hyperparameters used in Section 4.4 are listed for SWIM (Bolager al. 2023 Driver and rcond (12 reg.) are the parameters of t torch torch. www.llinnnnaa.... (Paskke et al. 201
Optimizer
Driver
Parameter sampler
Resample duplicates
1 reg.
SWIM Boluger et al.
2023
GELS
relu
Frue
1e-15
Table D.25: Hyperparameters used in Section 4.4 are listed.
Model
FGNN, FGNODE, LGN, HGN
GNODE, LGNN, HGNN
Node embedding dim.
Edge embedding dim.
hidden layers
2
hidden neurons (per layer)
16
message passing layers
1
Activation
squareplus Adam
squareplus Adam
Optimizer
Learning rate
3
10
Batch
100
Epochs
10000
10000
Precision
double
12 regularization
Table E.26:
Hardware used for the experiments is listed with details on CPUs (Intel i7 and AMD EPYC), memory, GPU (NVIDIA), CUDA version (driver version, CUDA version), and operating system (OS) versions of Ubuntu LTS, together with memory requirements (in GB).
Experiment
CPU (cores)
Memory
GPU (yrm)
CUDA
OS
section
17-14700K
66
RTX 4070 (12)
$50.120, 124
24.04.2
fig. fig. 4x4 train
rain
17-14700K
RTX (12)
12.4
17-14700K (20)
66
24.04.2
fig.
EPYC 7402 (24)
256
20.04.2
fig.
16
EPYC 7402 (24)
256
20.04.2
section
molecular systems
17-14700K (20)
66
24.04.2
section
EPYC 7402 (24)
256
20.04.2
Table F.27: Ablation study showing the influence of widths of the dense layers encoding the node and edge features and the linear layer on the mean squared error between predicted and true system dynamics y.
dg 128
dz = 256
dg = 512
da
2.879c-02
1.336c-02
3.666c-03
dg 16
16657-00
1.634e-04
6.571e-05
32
8,290e 04
2.011:004
1.278e
d, 64
7.826e-04
3.037e-05
3.713e-06
128
10
256
10 2
di 512
10
10
MSE
MSE
d
10-4
10-4
d, 16
dn 32
10-5
10-5
d, = 64
8
16
32
64
128
256 Linear layer width
512
Dense layer width
Figure F.7: Ablation study for widths of the dense and linear layers.
26

---

Published as a conference paper at ICLR 2026
Table F.28: Ablation study showing the condition number of the matrix associated with Equation ( with increasing feature widths.
128
dr. 256
dr. 512
4.266e+05
8.985e+06
2.612e+08
dg 16
1.525e+08
1.687ex06
3,545c407
dn 32
8.105e+04
1.352c+06
2.974c407
64
1.072ex05
9.625ex05
1.766e407
Table F.29:
Ablation study showing the influence of applying recursive message-passing on the mean squared error between predicted and true system dynamics ÿ. #msg is the number of message passes.
+msg
kmag=
#mag=
#mag=
#mag= 5
#mag=
Summing
7.713c-06
1.866-055
1.271c-04
1.150<-03
7.620c-04
2.223x-03
Averuging
2.020e-02
1.265e:02
1,280e-02
1.057e-02
1.175e002
1.089e002
as the system size (network width d,,) increases, but increasing the feature width of the encoders dN slightly stabilizes the system for large network widths, expected.
We additionally vary the number of message passes (#msg) by recursively applying the local-pooling hy ΣieN, NY aggregating the node encodings hy of the neighboring nodes i NN (#msg
times), and then applying the final message scheme explained in Section 3.2.2 with two different schemes: summing (hy + ΣiEN, h2 and averaging (hy + NN LIEN, h Table F.29 reveals
improve the accuracy for the 8- partttllee mass-spring system when averaging is used. We believe that summing works better than averaging because it implicitly encodes the node degree information by aggregating the neighboring messages. Each neighboring message is the output of a softplus activation function and has non-negative values. In all the other experiments presented in this paper, we only single message pass and do not optimize the number of message passes, as all the ground truth systems we consider only require single step of neighborhood information.
that having multiple message-passer
G COMPARISON WITH BENCHMARK DATASET
To further support our claims, here we perform benchmarking of our model against existing suitable graph network approaches. We made use of the existing publication from the NeurIPS 2022 Datasets and Benchmarks Track by (Thangamuthu considered models for comparison include:
(2022) and their corresponding repository. The
Full Graph Neural Network (FGNN) Based on the work of Saacaaez--nzzalee al. (2020), these models utilize message-passing as a key feature to enable simulation frame- work. Note that in the original work the architecture is called Graph Network-based Simulators (GNS) but for benchmarking it is called FGNN and we use this name as well.
Full Graph Neural ODE (FGNODE). An ODE version of FGNN is what we refer to as FGNODE (Sanchez-Gonzalez et al. 2019).
Graph Neural ODE (GNODE) This architecture uses graph topology to parameterize the force of a system using neural ODE approach, it was introduced by Thangamuthu et al.) (2022).
Lagrangian Graph Network (LGN) This architecture uses an FGNN to predict the Lagrangian of the system (Bhattoo et al. 2022).
Lagrangian Graph Neural Network (LGNN) Similar to LGN, this architecture decoples ( (Thangamuthu et al. 2022).
the kinetic and potential energies
Hamiltonian Graph Network (HGN) In this architecture FGNN predicts the Hamilto-
nian of the system (Sanchez-Gonzalez et al. 2019 Thangamuthu et al. (002))
Hamiltonian Graph Neural Netwrok (HGNN) Analogously, this architecture is similar to HGN but decouples the potential and kinetic energies of the Hamiltonian (Thangamuthu
et al.
2022
27

---

Published as a conference paper at ICLR 2026
Table G.30: Comparison the SOTA physics-informed graph network architectures (also see Table and our (SWIM) RF-HGN.
Model
(SWMM)
RF-
FGNN
FGNODE
GNODE
LGN
LGNN
HGN
HGNN
HGN
Translation invariance
X
X
X
X
X
X
Rotation invariance
X
X
X
X
X
Energy conservation
X
Gradient-deseen-ffice training
X
X
X
X
X
X
X
2x2 3D-Lattice System Loss
3x3 3D-Lattice System Loss
4x4 3D-Lattice System Loss
10a
[Adam] HCN
(Adami HGN train
[Adam) HGH train
107
RF-HGN
10r
[Adam)] last
(Adam) HGN bust
(Acam) HCN
MSE
(SWIM) RF-HGN best
MSE
107
(SWIM) RF-HGN test
MSE
(SWIH) RF-HGN test
10-2
10-4
10
10
0
2000
4000
6000
8000
10000
2000 4000
6000
8000
10000
2000
4000
6000
8000
10000
Training step
Training step
Training Step
Figure H.8: MSE losses the training and test dataset for a 2x2 (left), 3x3 (middle) and 4x4 (right) lattice during iterative training are given with solid lines for the average over ten runs; the shaded region extends from the minimum the maximum value. The dashed lines denote the (constant) MSE losses for our non-iterative optimization, and shaded dashes show the minimum and maximum.
First, we highlight the similarities and differences in model properties in table G.30 noting that our model satisfies requirements necessary for modeling physical systems while maintaining energy conservation.
A key metric of interest for our work is the training time, thus we have re-trained models from (Thangamuthu al. 2022) on their datasets of the spring system with 3, and nodes and record the training time. All runs were performed on the same machine as the experiments in the main paper. The resulting training times in table Table show clearly that our proposed approach much faster to train, especially compared to the specialized Lagrangian and Hamiltonian graph networks.
Of course, model is only useful if it can accurately make predictions, thus we plot errors on a test trajectory for all mentioned models in Figure H.22 We observe that (SWIM) RF-HGN has similar predictive ability as the SOTA architectures. It should be noted that for the test trajectory shown for N = 4 the LGN model diverged after around 25 steps. Similar results of diverging models were also observed from LGNN and the NODE architectures when we attempted to test on 100 trajectories, where multiple predicted trajectories would diverge from the true trajectory.
H ADDITIONAL RESULTS
Our submitted folder contains an animation of the test system shown in Figure of the main text, as well as the molecular dynamics systems (see Figure H.I8land Figure H.20
H.1 BENCHMARKING AGAINST SOTA OPTIMIZERS
In Figure H.8 we show loss curves for the Adam optimizer, highlighting how its train and test losses evolve over time relative to the loss of our non-iterative approach in 2x2, 3x3, and 4x4 systems. The model and optimizer hyperparameters are set accordingly as explained Appendix C.I Table D.7 Table D.9 and Table D.8 We observe comparable accuracies of Adam and SWIM (Bolager et al. 2023), even after 10000 gradient descent iterations using the Adam Kinggma & Ba 2015) optimizer. Moreover, Figure reveals that method scales better than iterative optimization, maintaining
low error as system size increases. And Table revells two to three orders of magnitude quicker training of (SWIM) RF-HGN than (Adam) HGN in different 3D lattice systems.
28

---

Published as a conference paper at ICLR 2026
(Adam)
10a
(Adam) HGN test
(SWIM) RF-HGN test
[Adam] HGN train
10-1
10a
Training time (5)
10a
Helobiee
10a
10"2
10-1
2X2
3
4
3
4 4
Lattice siee
Lattice sire
Figure H.9: Relative error and training time are shown for different lattice sizes. Boxplots show the mean and error bars based on ten runs with different random seeds.
Table H.31: Summary of the training times of the experiment presented in Figure H.9 for (SWIM) RF-HGN and (Adam) HGN in seconds. For the systems of sizes 2 (GPU trained), 3 (GPU trained), and 4 x (CPU trained), we observe approximately three, two, and three orders of magnitude faster training, respectively.
System size
(SWIM) RE-IGG
(Adam) HIN
2 x
0.06 seconds
Ra 82.93 seconds
x
0. seconds
92.1 seconds
4 X
3.06 seconds
936.12 seconds
using true 0.00
(SWIM) RF-HGN
(ELM) RF-HGN
(Adam) HGN
1.22
MYNY
0.5
2.12
of node
0.50
12-1
2/05
1.00
12-1
3.00
16-1
2002
7502
10603
2500
5200
7500
10002
2500
5000
7500
10200
2030
5000 7030
10000
Time
Time stop
Time stop
Time stop
Figure 1010: Illustration of position trajectories of the middle node over time (also see Figure H.11 Models are trained with nodes and tested (here) with 10 nodes using the standard spring potential (see Appendix C.I
wing
(EIN] AF-IGG
[Mdem)
33
10
28
55 50
Figure H.11: Illustration of position trajectories over time from models trained on system with nodes using the spring chain potential (see Appendix C.I and zero-shot tested with nodes with an external force (gravitational)
29

---

Published as a conference paper at ICLR 2026
using true X
(SWIM) BE HGN
(ELM) RF HGN
[Adam] HGN
14
7.20
-9.45
2.2
9.30
10
0.75
True
MSE
30
0.8
11.00
10"
10
-1.25
10
68
+1.50
2020
5020
7500
19532
2020
5000
7500
10000
2500
2000
7500
10000
2500
5000
7530
10000
Time stop
Time stop
Time stop
Time Stop
Figure H.12: Illustration of position trajectories the middle node over time (also see Figure 1. Models are trained with 5 nodes and tested (here) with 10 nodes using anharmonic spring potential (see Appendix C.I
wang
ISW
0020000024
C000
000000009
2. 85
3.0 1.9
200
2.
80 00
24
Figure H.13: Illustration of position trajectories over time from models trained on system with nodes using anharmonic spring potential (see Appendix C.T and zero-shot tested with 10 nodes with an external force (gravitational).
using true N
(SWIM) RE-HGN
(BILM) RP HG-HGN
(Adam) HGN
1.25
10-1
3.22
10
3.15
10m1
10
1:18
3.08
2300
5002
7300
12600
2000
2000
7502
20020
2000
5000
7500
20002
2500
5000
7500
20000
Time StOp
Time stop
Time stop
Time Stop
Figure 14 Illustration of position trajectories the middle node over time (also see Figure H.11 Models trained with 5 nodes and tested (here) with 10 nodes using the Morse potential (see Appendix C.T
000000000
-16-
ReoouooR
45.
48
16-
-15-
00
05
20
26
35
26
155
2.0
60
65
10
15
26
00
35
26
15
2. 65
85 20 1.5 20
Figure H. Illustratoon of position trajectories over time from models trained on a system with nodes using the Morse spring potential (see Appendix C.Ip and zero-shot tested with 10 nodes with an external force (gravitational).
30

---

Published as a conference paper at ICLR 2026
(MUM) HP-HIN
ISWINO Ri-MAN
(H RI-MUN
RF-MUN
Tainina time 14
10-1
10
100
10-1
10-4
210
212
Number of nodes
Number nodes
Figure H.16: Zero-shot generalization in 2D open chain (see Figure 4 (b)): RF-HGN trained up to N accurately generalizes up N 4096, outperforming retrained RF-HNN (right). RF-HGN with zero-shot generalization is also faster than RF-HNN for node counts larger than 2 (left).
HNN
(SIMM) RE-HNN
(SWIM)
4: of node 0
Qx of node 0
0.2
6.2
0.1
of node
8.0
rode
0.0
-0.2
-02-
4000
2003
4000
3000
4000
2000
1000
step
Time
Time
Time step
0.33
0.4
0.25
9t of nade
node
Ot of
0.0
0.00
-04
-025
0.00
2002
8200
2003
8000
2000
4000
2000
4003
Time step
Time step
Time step
Time step
Figure H.17: Illustration of position trajectories of the corner nodes over time (also see Figure 6 Top: Trained with 2 nodes and tested with 2 nodes. Bottom: Trained with 2 nodes and tested with 212 nodes.
H.2 ZERO-SHOT GENERALIZATION AND COMPARISON OF RANDOM FEATURE METHODS
Figure H.17
illustrates the trajectories of the corner nodes of the experiment in Figure For this particular example the left corner trajectory (node with id 0) seems to be harder capture than the other nodes in the system for the extreme zero-shot case (trained with 8 nodes, tested with 4096 nodes) case, hence slightly higher error the trajectories compared to the non-zero-shot case (trained with nodes and tested with 8 nodes) as one can see in Figure
H.3 BENCHMARKING AGAINST SOTA ARCHITECTURES
H.4 ROBUSTNESS AGAINST NOISE
To evaluate model's robustness against additive noise, add Gaussian noise with different standard deviations (σ) to the ground truth positions and momenta before including them in the training set. We experimented with nodes in 3D the training set with 1000 samples in phase
Table H.32: Molecular dynamics evaluation with 9 particles. Mean squared error (MSE) and relative 12 error (rel. (2) are reported together with the true Hamiltonian over the ground-truth trajectory and the (ELM) RF-HGN predicted quantity over the rolled-out trajectory.
T=1
T=25000
T=50000
T=74999
T=99999
q rel. 12
8.651e-08
9. 9.369e+07
2.926e+09
9.456er08
6.337ex09
2.044c-04
6.766e+03
3.718c+04
2.139s404
5.615e404
Thue H
-1.333+01
-1.233e=01
-1,233e+01
-1,233e+01
.1.233e+01
Modal R
-14450e+01
2.979e+10
-1.122e+10
-3.085e+10
8.765e+10
31

---

Published as a conference paper at ICLR 2026
Table H.33: Molecular dynamics evaluation with 9 particles. Mean squared error (MSE) and relative 12 error (rel. (2) are reported together with the Hamiltonian over the ground-trute trajectory and the (Adam) HGN predicted quantity over the rolled-out trajectory.
T=1
T=25000
T=50000
T=74999
T=99999
q
5.043e.03
1.554e-02
6.052e-02
1.279-001
rel. 12
3.863e-07
4.94e-02
8,568c-02
1.711e-01
2.522c-01
True
-1.233c401
-1.233c=+11
-1,233c+01
-1.233c+01
-1.233c+01
Model R
-1.233e401
-1,233c=01
-1,233e+01
-1,233e+01
-1.233c+01
26-
25-
Treestes
Time step 2500e
50893
70599
59999
Figure H.18: Illustration of position trajectories time from models trained system with nodes on the 2D molecular dynamics system (see F (c)).
using
1.4
dr of node
True Kt
20-1
12-
10-1
of
1.2
10-1
08-
1.0
10ite secto 15002 305007
23100 stoos roooe 200302 Time step
25200 suo2e rstoo 300500
25030 52tce 15t00 120000 Time step
Time
Time
Figure 19: Illustration of position trajectories over time from models trained on a system with nodes the 2D molecular dynamics system (see F (c)). Results from ELM RF-HGN training are omitted due to very large errors which distort the representations the plots.
32

---

Published as a conference paper at ICLR 2026
(Adem)-IGU
RF-HIG
Tiu 25000
Tile 50030
Tree 3000s
Figure H.20: Illustration of position trajectories time from models trained system with 36 nodes and zero-shot tested with nodes on the 2D molecular dynamics system (see Figure (d)).
0.6
using true N
(5W/M) RF-HGN
[Adrm)] HGN
1.7
0.4
10
131
of node 32
-24
42-
10
-36-
40
-101
10
10-1
200003
130002
step
Time step
stop
Time step
Figure H.21: Illustration of position trajectories time from models trained system with 36 nodes and tested with 64 nodes on the 2D molecular dynamics system (see Figure (c)). Results from ELM RF-HGN training omitted due very large errors which distort the representation in the plots.
RF-HGN
FGNN
PENODE
GNODE
HGN
HGNN
LGN
LGNN
N=3
N=4
N=5
10-1
Bollout Entor
10-1
10-2
10
10
20
40
80
100
20
80
100
40
100
Time Step
Time Step
Time Step
Figure H.22: Rollout errors on test trajectory benchmark for N 3, 4, 5.
33

---

Published as a conference paper at ICLR 2026
Table H.34: Results of training with noisy data are displayed. Relative 12 and mean-squared errors are displayed.
le:5
1e-2
CV 6.18c-3
6.59c-3
5,35o-3
7.84o-
6.45c-3
3.38o-2
1.1e+1
CV 12
4.03e-2
4.07c-2
3.72c-2
4.29-2
4.05e-2
9.61e-2
6.64e-1
8.78e-1
Test MSE
9.52e-3
9.77e:
9.12e:3
1.0e-2
9.12e:3
9.67e:3
1.01e-1
3.94e-1
Test
4,96e
4.92e-2
4.83e-2
5.09e
4.80e-2
4.94
1.65e
3.27e-1
Table H.35: Batch-wise training results compared direct least squares solutions displayed with training time in seconds and memory usage in GiB. Relative 12 and mean-squared errors are displayed.
ELM
ELM (batched)
SWIM
SWIM (balched)
CV MSE
9.01
2. 1.e5+003
3.84e-1
1.08
CV rel.
2.59c-1
2.54
4.65e-2
7.89c-2
Test MSE
1.98+33
3.71c-1
1.09
Test rel. 1
2.60-1
2.49
4.61e-2
8.02e-2
Trinnin time
3.96
4,48
3.88
4.42
Memory usage
2.4
1.1
2.4
1.1
space using 5-fold cross-validation, and on test set of 1000 samples in phase space. To stabilize the results, we repeated this experiment times with different seeds and report average results in Table H.34
The columns in Table H..4 list the standard deviation a used the experiments. The results indicate the robustness of our method to Gaussian noise added to the state. For more realistic scenarios, the model can further be improved by modeling uncertainties and assessing sensitivities for noise in real-data setting.
H.S BATCH-WISE WISE TRAINING
To demonstrate batching, we prepared an example with Lennard-Jones potential with 4 particles in 2D trained with both ELM, SWIM, and their batch-wise versions when solving the linear problem (Equation (6)) Batching is performed by sub-sampling the training data set (5000 states in phase space) and averaging the resulting last-layer coefficients. Table H.35 lists the results, where the columns indicate which random feature method is used, rows indicate cross volsd-aaliatiin (5-fold) errors with train set size 5000 and test errors averaged over different seeds of test size 5000.
We note that the memory required for batch-wise training can be further improved with additional tuning, for example, by explicitly freeing the GPU memory, we could tune the memory requirement to be as minimal as possible (0.04 GiB). In this case, however, the training time also increases to around 5.7 seconds. Different tuning strategies (e.g., for lower memory, for quicker runtime) are therefore important consider when comparing different training strategies.
We believe that with an established linear solver like LSQR (Paige & Saunders 1982. LSMR (Fong & Saunders 2011), LSRN (Meng al.) [ (2114)), one can further study the batch-wiss training of our model for even larger systems, tuned for specific needs such as low memory or fast training.
H.6 BENCHMARKING DIFFERENT RANDOM FEATURES
We also experimented with random Fourier features (RFF) (Rahimi & Recht 2007) by setting
Wi ~ NNO,, ) ~ Uni Unifform(( 2m), z #feeturrs cos Cos(WT., Tx x + b),
where: the random features and features the size of
We have extended the noise-scale experiment in Appendix H.4 with R = 1 and list the results in table Table H.36 Additionally, we run the same experiment using ELM by setting
Wi ~ (' 1), b Uni Siniform(-T,r
and list the results in Table H.37
34

---

Published as a conference paper at ICLR 2026
Table H.36: Results of training with noisy data using RFF are displayed. squared enrors displayed.
Relative 12 and mean-
1o-4
le-1
MSE
3.03
3.05
3.01
3.09
1.4++1
CV rel. (2
9.1e-1
9.14e-1
9.15c-1
9.13e-1
9.1e-1
9.17c-1
9.41e1
9.87e-1
Test MSE
3.43
3.61
3.44
3.47
3.48
3,39
3.62
3.87
Test rel.
9.74e
9.77e-1
9.82e-1
9.84e
9.71e-1
1.04
Table H.37: Results of training with noisy data using ELM are displayed. Relative 1 and mean- squared enrors displayed.
le.5
le:3
le-2
CV MSE
3e-1
3.2
2.44c-1
6.78e-1
2.4c-1
4.62e-1
8.75
2.32e+2
f2
2.54e-1
5.23e-1
2.39c-1
3.24e-1
2.31e-1
2.87e-1
8:73e-1
Test MSE
4.7c-1
2.89c-1
2.63c-1
3.65e-1
3.04c-1
3.06e-1
5.91++1
2,42046
Test rel. 1
3.08e-1
2.76e-1
2.66c-1
3.02e-1
2.76e-1
2.76c-1
1.4
1.65+42
We note that data-agnosia methods perform better when tuned slightly towards the problem. To demonstrate this, we experimented on the same system but without noise. Table H.388) we list relative 12 errors for RFF; the parameters ofF and oF are node/edge encoder and message encoder RFF parameters, respectively. We note better approximations with larger sigmas (lower standard deviation), similar to what we have observed with ELM. This is an important point, highlighting the value data-agnostic methods in certain cases, especially when their additional tunable parameters set appropriately. In our main paper experiments, mainly chose the best-performing random feature method that requires no extra tuning, enabling fast training while maintaining accuracy comparable to tradient-descent- based approaches. However, in real-worke scenarios, this fast training could be leveraged to further tune the hyperparameters of RFF or ELM and select the best-performing configuration.
Table H.38: RFF results are displayed where of is the node and edge encoder RFF parameter, and of is the message encoder RFF parameter. Relative 12 errors are displayed.
02
0.1
02
0.1
1.02
1.01
1.01
6.52c-1
10
5.34e-
7.04c-2
4,76e-2
35


Reviews:
[{"review_id": "pUUhyJqhkf", "reviewer": "Reviewer_RhM6", "summary": "This paper proposes RF-HGN, a method to accelerate training of Hamiltonian Graph Networks by replacing gradient-based optimization with random feature sampling and least-squares solvers. The approach samples dense layer parameters randomly (using ELM or SWIM) and optimizes only the final linear layer, achieving 100-600× speedups while maintaining comparable accuracy. The method demonstrates zero-shot generalization from small training systems (8 nodes) to large test systems (4096 nodes) across mass-spring, lattice, and molecular dynamics systems.", "questions": "1. Can you provide convergence guarantees or approximation bounds for the random feature approach in the\n  physics-informed setting?\n2. Why specifically Hamiltonian neural networks? How does performance compare when applying random features to\n  Lagrangian or other physics-informed architectures?\n3.Can you provide principled guidelines for when accuracy degradation becomes significant? What physical properties are most affected?\n4. What are the practical limits of your approach? At what system size does the linear solver become prohibitive?\n5. How does the method perform on more challenging physical systems beyond simple spring-mass dynamics?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. Dramatic Practical Speedups: The 100-600× training acceleration is genuinely impressive and could enable new applications in physics\n  simulation.\n2. Strong Zero-Shot Generalization: Training on 8-node systems and testing on 4096-node systems demonstrates remarkable scalability - this is perhaps the paper's most valuable contribution.\n3. Comprehensive Experimental Validation:\n    - Comparison against 15 different optimizers provides robust baselines\n    - Multiple physical systems (springs, lattices, molecular dynamics)\n    - Use of established NeurIPS 2022 benchmark dataset\n4. Physical Consistency: The method maintains energy conservation and incorporates essential symmetries (translation, rotation, permutation invariance).\n5. Clear Algorithmic Contribution: The two-stage training procedure (random sampling + least squares) is well-defined and reproducible.", "weaknesses": "1. Weak Theoretical Foundation:\n    - No convergence guarantees or approximation bounds\n    - Limited analysis of when/why the method works\n    - Missing connection to random feature theory for this specific setting\n2. Poorly Motivated Architectural Choices:\n    - No compelling justification for choosing HNNs over alternatives (Lagrangian, Port-Hamiltonian, etc.)\n    - Graph networks not well-motivated for many test systems (regular lattices better suited for CNNs)\n    - Missing literature review of Hamilton graph neural networks\n3. Limited System Complexity:\n    - Mostly simple spring-mass systems and basic molecular dynamics\n    - ~10% relative error on Lennard-Jones systems suggests limitations for complex potentials\n    - No testing on truly challenging physics (e.g., turbulence, phase transitions)\n4. Accuracy Trade-offs Not Well Characterized:\n    - Sometimes less accurate than second-order methods (LBFGS)\n    - Large variance in results (Table 1) raises questions about reliability\n    - No principled way to predict accuracy vs. speed trade-offs\n5. Scalability Questions:\n    - Memory complexity O(MNe) may become prohibitive for very large systems\n    - Linear solver bottleneck O(Kd²L) not thoroughly analyzed\n    - Integration constant handling appears ad-hoc\n6. Missing Key Comparisons:\n    - No comparison with other physics-informed ML acceleration techniques\n    - No baseline comparisons with structure-specific alternatives (CNNs for lattices)\n    - Limited comparison with other random feature applications to physics", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "uivBRpEohA", "reviewer": "Reviewer_5YK6", "summary": "This paper introduces RF-HGN, a method for training Hamiltonian Graph Networks using random feature sampling instead of iterative gradient-descent optimization. The authors demonstrate significant training speedups compared to other optimizers (Adam, LBFGS, etc.) while maintaining competitive accuracy on mass-spring and molecular dynamics systems.", "questions": "- The method is presented in the context of Hamiltonian systems. How readily can it be applied to other physics-informed graph networks, such as Lagrangian or Port-Hamiltonian networks, or even non-conservative systems? Does the \"random features + linear solve\" recipe generalize? \n\n-  The method demonstrates impressive zero-shot generalization to larger systems, but does this generalization hold for structurally heterogeneous systems? For instance, if a model is trained on a regular lattice (where all nodes have the same degree), can it accurately predict the dynamics for a system with a mix of node degrees, or for a node with an unexpectedly high degree not seen during training?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "The proposed method offers a significant speedup without sacrificing accuracy in multiple systems. \n\n- comprehensive benchmarking\n- elegant solution by merging  Hamiltonian Graph Networks (for physics-informed modeling), Random Features (for fast, non-iterative training), and careful construction of physical invariances (translation, rotation, permutation). \n- generalizability from small to larger systems", "weaknesses": "- the paper could give more background on why random features work\n- acknowledge the limitation more prominently and discuss it\n- the acceleration is significant, but the claimed 600× speed-up is overstated; compared to the second-best model, it is actually 150×.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IxDurxACAQ", "reviewer": "Reviewer_258K", "summary": "The paper introduces Random-Feature Hamiltonian Graph Networks, whose hidden layers are fixed random features, with only the final linear readout solved in one least-squares step instead of iterative gradient descent. The model encodes translation/rotation/permutation invariances, reports 100-600x training time speedups over many optimizers on mass-spring and Lennard-Jones systems, and shows zero-shot scaling from tiny to very large graphs without retraining.", "questions": "1. What are the key differences between your approach and reservoir computing? \n2. What happens with 2-3 random feature message passing blocks versus a single block? How does performance scale with feature width, and do you observe conditioning issues in the least-squares solve as capacity grows?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. Replacing long iterative GD with a single convex least-squares solve is simple and yields substantial speedups without complicated tuning.\n2. The model encodes translation, rotation, and permutation invariances appropriate for N-body dynamics, which improves data efficiency and generalization within a family.\n3. Demonstrates training on small systems and inference on much larger ones without retraining.", "weaknesses": "1. The core idea is somehow similar to reservoir computing: both avoid training the feature generation module and optimize only the final readout layer.\n2. With a single message passing stage and a linear head, it’s unclear how the method scales to long range or multi scale interactions. There is no ablation on stacking multiple RF blocks.\n3. The training formulation assumes no external forces and exact energy conservation. Robustness to mild non-conservation (damping, stochasticity) is only lightly tested.\n4. Transfer across distinct topologies is under-explored.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "S4kwSDgPWd", "reviewer": "Reviewer_nwfX", "summary": "This paper proposes a novel iterative method to train Hamiltonian Graph Neural Networks (GNNs) orders of magnitude faster than standard optimizers. The approach leverages random feature sampling over the small shared node, edge, and message MLPs of the GNN, followed by a least-squares fit to determine the final linear layer that outputs the Hamiltonian. To ensure translation and rotation invariance, the particle positions are processed by aligning the simulation frame of reference to the center of mass and rotating the coordinates to a fixed orthonormal basis. The proposed training algorithm is evaluated on four different graph configurations and compared against multiple standard optimizers, achieving at least two orders of magnitude speedup (even outperforming the second-order L-BFGS optimizer). Furthermore, the method is benchmarked against seven state-of-the-art structure-preserving architectures using publicly available datasets.", "questions": "* Section 4.3: Is there a particular reason why coosing Hamiltonian GNNs instead of, say, Lagrangian GNNs? Have the authors tried to train the tested architectures (GNODE, LGNN, FGNN, etc) with the presented method to see which architecture performs better?\n* Have the authors tried to use simpler invariance tricks, like using relative distances or relative angles?\n* Line 1079: $ 5^{-3}$ might refer to $5\\cdot 10^{-3}$?\n* Line 1079: Here is specified that the molecular dynamics experiments are only rolled-out to 50 timesteps. However, Table 2 shows T=99999 timesteps. Which is the correct one? 50 timesteps is very few time horizon for a molecular system, given that the timestep is very small for stability.\n* Equation 6: Is the least squares minimization well-posed? Have the authors found any problems when training? I'm thinking about a bad conditioning number for the normal equations $Z^TZ$, or some degeneracies induced by the rotation-translation invariance.", "limitations": "", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "* The method is tested over Neurips 2022 open source benchmarks.\n* The paper proposes an interesting way to enforce rotation-invariances in the system.\n* The experiments are varied, exploring N-body systems (chains, regular grids) and molecular interactions in Lennard-Jones systems\n* The generalization and rollout results are good. These are the precise experiments needed to show the power of structure-preservation in learning conservative dynamics.", "weaknesses": "* Even though the method is tested over a wide variety of systems, they are still small scale toy problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)