Hamiltonian Neural PDE Solvers through Functional Approximation
Authors: Anthony Zhou
Amir Barati Farimani
Venue: NeurIPS 2025 poster ()
Keywords: Neural PDE Solvers
Functionals
Neural Operators
Hamiltonian Neural Networks
Primary area: machine_learning_for_sciences
TL;DR: A new method for predicting PDE solutions is proposed based on Hamiltonian mechanics and functional approximation.

Abstract: Designing neural networks within a Hamiltonian framework offers a principled way to ensure that conservation laws are respected in physical systems. While promising, these capabilities have been largely limited to discrete, analytically solvable systems. In contrast, many physical phenomena are governed by PDEs, which govern infinite-dimensional fields through Hamiltonian functionals and their functional derivatives. Building on prior work, we represent the Hamiltonian functional as a kernel integral parameterized by a neural field, enabling learnable function-to-scalar mappings and the use of automatic differentiation to calculate functional derivatives. This allows for an extension of Hamiltonian mechanics to neural PDE solvers by predicting a functional and learning in the gradient domain. We show that the resulting Hamiltonian Neural Solver (HNS) can be an effective surrogate model through improved stability and conserving energy-like quantities across 1D and 2D PDEs. This ability to respect conservation laws also allows HNS models to better generalize to longer time horizons or unseen initial conditions.

Hamiltonian Neural PDE Solvers through Functional Approximation
Anthony Zhou
Amir Barati Farimani Carnegie Mellon University
Carnegie Mellon University ayz20andrew. cmu. edu
barati@cmu.ed edu
Abstract
Designing neural networks within Hamiltonian framework offers a principled way to ensure that conservation laws are respected in physical systems. While promising, these capabilities have been largely limited to discrete, analytically solvable systems. In contrast, many physical phenomena are governed by PDEs, which govern infiniteddimensionall fields through Hamiltonian functionals and their functional derivatives. Building on prior work, we represent the Hamiltonian functional as a kernel integral parameterized by a neural field, enabling learnable sunction-to-scalle mappings and the use of automatic differentiation to calculate functional derivatives. This allows for an extension of Hamiltonian mechanics to neural PDE solvers by predicting functional and learning in the gradient domain. We show that the resulting Hamiltonian Neural Solver (HNS) can be an effective surrogate model through improved stability and conserving energy-like quantities across 1D and 2D PDEs. This ability to respect conservation laws also allows HNS models to better generalize to longer time horizons unseen initial conditions.
Introduction
Physics is remarkably complex; a of governing laws can produce endlessly diverse and interesting phenomena. In effort unify governing laws, Hamiltonian mechanics describes physical phe- nomena under the Principle Least Action, whereby the trajectory system evolves because has locally the least action. Remarkably, the minimization of action can derive nearly all physical phenomena, from kinematics to quantum mechanics, suggesting an underlying, fundamental property of nature. Rather than describing physics through forces and masses, such as in Newtonian mechanics, using energy and symmetries allows for an elegant and unifying perspective of physical systems.
Despite the long history of Hamiltonian mechanics, most current numerical and neural PDE solvers operate within Newtonian framework. This has been very successful; forces and masses are easy to intuitively understand, and both numerical and neural solvers have found important, practical uses. However, motivated by natural conservation laws and symmetries that arise from Hamiltonian mechanics, we seek to investigate neural PDE solvers that operate in this framework. We find that the Hamiltonian framework brings inductive biases that improve model performance, especially when conserving energy-like quantities.
Extending Hamiltonian mechanics neural PDE solvers requires both theoretical and implementation insights. In particular, most PDEs are used to describe continuum systems (fluids, waves, elastic bodies, etc.); as such, infinite-dimensiona Hamiltonian mechanics are needed. this setup, the Hamiltonian is defined as functional that maps from input fields to a scalar, often interpreted as the energy of the system. The variational or functional derivative of the Hamiltonian then contains dynamical information that evolves the system. This introduces interesting modeling challenges; in particular, the need approximate infinitediimensionll scalar mappings and for the approximation to have functional derivatives.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In this paper, we consider models based on learnable kernel integral, framework that can provably approximate linear functionals and implicitly learn functional derivatives. Through automatic dif- ferentiation, functional derivatives can be obtained and optimized in Hamiltonian framework to predict future PDE states. This allows the development of Hamiltonian Neural Solvers (HNS), which we evaluate on three PDE systems and find both high accuracy and generalization capabilities, which we hypothesize arise from inductive biases to conserve energy. Code and datasets for this work are released at https: //tp:th/ggttttt// com/ anthonyzhouu http:/hhhooohaammtttiiiieoeee _o_
2 Background
Hamiltonian Mechanics Hamiltonian mechanics is usually introduced in the discrete setting, where there are set of n bodies, each with a position and momentum. Therefore, the state of the system can be described by quantities, usually assembled into position and momentum vector P.P E R". In the discrete case, the Hamiltonian (H(q, p) 2 RR is function that takes two vectors and returns a scalar H. usually interpreted as the energy or another conserved quantity. Evolving the system in time is done through Hamilton's 's equations, often expressed in terms of the symplectic matrix and state vector u.
dq dt
OH 0p"
dp
OH
du /d. = dp/dt] = 0 T
8H/88q 0H/OP
= JV.H
(1)
dt
dq
dt
Hamilton's Equations
Symplectic Form
Derivatives can be evaluated from vector calculus identities; for clarity we consider the state vector with 2n components and write du ddul duan dt
= dt
and V.? Du
8H
8H
de
Juz.
Dnfinite-Dimensionall Systems Many PDEs describe continuum systems, such as waves, fluids, and elastic bodies, where the positions and momenta of discrete bodies not well defined. Despite this change, continuum systems can still conserve energy and be viewed a Hamiltonian framework. Firstly, finite-dimeniionl vectors q. p T become generalized functions u(æ, t) defined continuously varying coordinates x E Ω and at time The function u E F(2)) describes the state of the system by assigning a scalar vector for each coordinate (such velocity, pressure, etc.). Secondly, the Hamiltonian is generalized from a function (H(q. p) 2 + R} to functional [H(u]: : F(2) - RR} that returns scalar given input function. Lastly, the symplectic matrix generalizes to a linear operator J, and the gradient Vn becomes the variational or functional derivative 5 e F(2). This leads to Hamilton's is equations for dinntt--dmeesional systems, where F(2) is arbitrary test function [T]
du
Ou = 0t
8H1 du
8H
H(u he] H [1]
(z)e(z)dz= lim du
h-00 h
Hamilton's is Equations
Functional Derivative
From this background we can see why neural networks are natural choice for modeling discrete Hamiltonian systems: in this setting, the Hamiltonian function and its gradient naturally calculated with automatic differentiation. In the infinit-ddimensional setting, learning the Hamiltonian functional from data becomes less obvious.
Functional Approximation
Theory The development of functional approximators 2 can be seen a result of the Riesz representation theorem, a fundamental result functional analysis (6.. We give some intuition here, but provide full proof of linear functional approximation in Appendix [C.] We start by observing that a functional can be equivalent the inner product over suitable functions:
Theorem 3.1 (Riesz representation theorem). Let H be Hilbert space whose inner product (ar, y) defined. For every continuous linear functional up E H* there exists unique function J E H, called the Riesz representation of such that:
xpズ] (z)
for all H.
(3)

---

Immediately, we can see that the Hamiltonian M(u can potentially be expressed as an inner product Kg) over the imput function and learnable function Kig. This changes the problem of approxi- mating a functional to the problem of approximating a function, which neural networks can readily achieve through universal approximation theorems 7 g
A common inner product for arbitrary functions u, in a Hilbert space H involves an integral: (u,v) Jo u(z)v(z)dr. One can check that this satisfies the basic properties of inner products and is also valid for vector-valued functions V. Based on this formulation, functionals can be parameterized by an Integral Kernel Functional (IKF):
Holu] (o(x)u(()dzz
In this form, we observe connection between the integral kernel functional and the integral kernel operator that neural operators are built on [9]. Specifically, the integral kernel operator is given by:
(a(u)(x)= / (e y)u(y)dy
(5)
where a learnable operator K. acting on an input function u(x) E U2 produces output function ((()) (a)) E V for Banach spaces U, V: this expression is equivalent to the IKF when evaluated at a single point z. We can make similar modification to neural operators where the kernel is allowed to be nonlinear by giving function values as input: 60(x, y, u(r), u(y)). In our setting, we can adopt a nonlinear kernel by allowing the kernel to change based on the input function: (o, u(x)).
One shortcoming of this framework is that the approximation nonlinear functionals not proven. Most Hamiltonians are nonlinear with respect to however, we empirically find that the proposed architecture can still approximate nonlinear functionals by using a nonlinear kernel Ki (zr, u(x)). There are some potential hypotheses for this. The Riesz representation theorem is strong statement, implying that for linear functionals, there is a single, unique Ki E H that can be used with any u E H to approximate H[u]]. In deep learning we are less constrained; we can potentially allow (a(u) to depend on u and even relax the uniqueness of Ke to still be able to approximate H H[u]] Furthermore, under assumptions of regularity and compact support, for arbitrary functionals H we can always show that a function f H exists such that H[u] (u, f) by constructing f = [|
M|u]
u, with the
norm: ||ul|2 = / SS u(x)dde.
Implementation
There are two main implementation choices: upproximating the integral and parameterizing the kernel. The integral can be approximated by a Riemann sum:
Ho [u] = = (z u(z))u(z)dz =
Σ (o, ((z:)))(
(6)
with quadrature weights μ4. Interestingly, the full Riemann sum can be computed since the sum is
only calculated once; in neural operator literature, approximating the integral with a full Riemann sum is usually too costly since the sum is evaluated for every query point. This results in approaches that truncate the sum [10 or represent in Fourier space through the convolution theorem [13]. In our implementation, to accurately estimate the integral.
evaluate different quadratures but find that the trapezoidal rule is enough
To parameterize the kernel, we take the perspective from the Riesz representation theorem. In the linear case, the kernel is simply a function k²(x) E H that maps input coordinates to values in the codomain. This is exactly what neural field is; therefore, we can inherit the substantial prior work on neural fields [4] to effectively parameterize Kg.
One useful concept is the conditional neural field. In particular, the nonlinear kernel (a u(x)) can be seen as a conditional neural field whereby the kernel output changes based on the input function u(x). This concept can also be further extended to distinguish between local and global conditioning [5]: in neural field literature, the field can be conditioned on local information u, = u(x;;) or global information [u(#o),, ... u(rr))). Conditioning mechanisms can also be implemented in various ways. In the simplest case, the conditional information u, or u is concatenated to the input coordinate [16H18) however, we adopt approach based on Feature-wind Linear Modulation (FiLM) [19]
3

---

IKF Approximation:
Kernel Parameterization: 2,)
sin (Wr+bb)
Scale/ /Sift
uo 14,,,440 1,
AX
xb
In U1
o(Wu+bb)
xL
7/00
,000
Figure 1: Left: The functional H(ux is approximated by an integral kernel functional (IKF). The integral is further approximated by a Riemann sum with a learnable kernel Kø. Right: The kernel can be parameterized with various architectures and conditioning mechanisms. A local, FiL.Mcconitiooeed SIREN kernel is shown, where each SIREN layer has a dedicated FiLM network to produce scale and shift parameters. The output 2; for single point 24, which needs to be calculated Vi (l..... n
Beyond conditioning mechanisms, an important consideration is the architecture of the kernel itself. powerful neural field architecture the sinusoidal representation network or SIREN [20]. SIRENs have shown good performance not only in fitting neural field Ke, but also in accurately representing gradients Vkg. which may help fitting the functional derivative 8H/6u. With these architectural
(1) with local conditioning can be written as: choices, layer of the kernel Ki
(1) Fig (Is, u4) 7o (u:))in (Waa z; + b) + A6" (us)
(t)
where 7 and 99 are the FiLM scale and shift networks at layer I. Each SIREN layer has its own FiLM network and each FiLM network can also have multiple layers. Inputs F2, u, can optionally be lifted and the output z; can be projected back to the function dimension. Given the hierarchical nature of the proposed implementation, we also illustrate the architecture in Figure
final consideration during implementation batching. The layer update Equation is written in pointwise form, but we can observe that the computation does not depend on other points {(z}, uy) j * i}. and the weights are shared between points. Therefore, the forward and backward pass for each point can be done in parallel; in practice, the spatial dimension is reshaped into the batch dimension for efficient training and inference.
Hamiltonian Neural Solvers
So far, integral kernel functionals are largely agnostic to the down- stream task, only requiring that the input function u be discretized at a set of points (z,, u,). To use functional approximators as PDE solver, introduce the Hamiltonian Neural Solver (HNS). The time evolution of PDEs can be described by Hamilton's 's equations, which depend on the operator I and the functional derivative 8H/ /8u. In Hamiltonian PDE systems, J is analytically derived or can be looked up, is guaranteed to be linear, and is usually simple. A common example in 1D is choose approximate J with a finite difference scheme, although it can
J O,. Therefore, also be learned with neural operator. Lastly, we find that 8H/55u can be computed directly with automatic differentiation. This is empirically shown in Section [.1 but one way to see this is by considering the linear IKF. The gradient V. Ei= Ko (a;z)); 14; Axr (x)))))..... (o( )], up to a constant factor, which the discretization of the learned function (z))). In this case, the functional derivative is approximated by an arbitrary function, which is the desired behavior.
Within this framework, the training is described in Algorithm| For Hamiltonian systems, we assume knowledge of the analytical form of H and 8'H/Su, and labels for 5'H 8u or du/de dt can
Algorithm 1 Training HNS
1: repeat
2: Ho Σ (o, uz)us 1r
be calculated from training data using numer-
3: 6Ha
t- autograd (Ho, u)
ou
ical schemes. The training loss can either be evaluated in the functional form 8HHo/ /du or the temporal form J (8HH// /ōu); we implement the former since it more directly optimizes Ki.
4: L C= _ 5H
or || áHa du
du
de
5: +- Update(#, VoL)
6: until converged
During inference, the current state u' and coordinates x are used in a forward pass to compute Ha. A backward pass is then used calculate 8Ha/Su; using the operator J, prediction for du J &He is calculated. The estimated derivative is used to update the state n'+1 ODFint(u's dt dt
1==
du using a numerical integrator. In our implementation, use a 2nd-order Adams- Bashforth method.

---

Base
OOD (c, E [1,3])
Disc.
(zy E [-2, 2)
Metric
MLP
FNO
IKF
MLP
FNO
IKF
MLP
FNO
IKF
Fi Fi[u] 2.47c-5
2.76e-4
3,00c-16
0.046
0.268
2.13e-7
49.3
49.7
0.737
/Fi/ôu 0.083
0,066
1.15e-3
0.114
0.122
1.15e-3
2.64
2.70
0.077
Fuluu] 0.029
0.016
2.05e-3
6709
6493
2908
381
322
55,4
FFu//uu 1.33
2.10
0.045
1730
1723
1079
92.1
88.1
35.5
Table 1: Validation MSE is reported across different experiments with linear and nonlinear functionals. Errors calculated in the scalar Flu] and gradient domains 8FF fu to evaluate performance in fitting functionals and implicitly learning functional derivatives. Additionally, generalization to OOD inputs and unseen discretizations evaluated. Parameter counts are: MLP (7.5K), FNO (10.9K) and IKF (4.3K); each experiment is repeated 6 times and errors are reported as the average across seeds.
4 Experiments
4.1 Toy Examples
To better understand and situate integral kernel functionals, we examine their ability to model simple, analytically constructed functionals and compare their performance to other architectures. The experimental setup as follows: given a functional Flu), random polynomials u(x) cg.
c1."-1+.. ..p......ppp1p1 1p- are generated by uniformly sampling {c₂ E [a, b) =0.... .P). Therefore, the dataset consists of N pairs of polynomials and evaluated functionals Fl", for ......, .N. In practice, each function v (xxx) is discretized at set of points {z; i M,,
such that it represented as u [u" (x1), ..., un wr (xL)], and Flu") remains a scalar.
We consider two cases: linear functional Fi[u] - u(x)) x and a nonlinear functional Ful[u] - S (u(z)))&d ddr. In both cases, constructed polynomials v²(x) can be substituted and the definite integral can be analytically solved to use as training labels. Additionally, the functional derivatives are 6.F: = x² and &FF.. 3u², which are used to evaluate the gradient performance of deg models. We restrict the degree of u" (z) p 2 and sample c; E [- 100 training and 10 validation samples are generated for each case. Coordinates T; are uniformly spaced on the domain
8u
[-1,]] with M 100 points.
Two baselines are considered: (1) a MLP that takes concatenated inputs [x, u] RRM and projects to output Fg R through hidden layers, and (2) FNO with mean-pooling at the final layer to evaluate an operator baseline. The FNO baseline takes inputs with spatial and channel dimensions, or stacked inputs [x, u| E IR x2 and projects to output field E RR which is then pooled to Fe e R. For simplicity, the integral kernel functional (IKF) is used with a local MLP kernel and is linear for Fi and nonlinear for Fut- In the nonlinear case, FiLM is not used and [x,, u.] IR2 is concatenated as input to Kg. Lastly, models are trained with the loss L = Em [F(u((a)) Fo x)] training is supervised with the functional to evaluate if models can implicitly learn accurate derivatives.
To further evaluate generalization, after training we consider two cases: (1) OOD, where validation in- functions, (2) Dise., where validation inputs u(x) are uniformly discretized on larger domain x, 21 2] with a larger Ar. while keeping M = 100. In each experiment, validation errors are reported with the metric NV >. [|F[u"] Fo(" x) or N >. |1 autograd( Fo, u²) to measure performance in fitting functionals and learning accurate functional derivatives.
puts u(z) are sampled with polynomial coefficients E [1,3] to model
Results are given in Table We observe that architectures based on an integral kernel can accurately represent functionals and their derivatives, outperforming other architectures. Despite only training on scalars F [[]], IKFs can implicitly learn the correct functional derivatives, as plotted in Figure This is new; even on toy problems, traditional architectures do not have well-behaved functional derivatives, as such, learning in the gradient domain for Hamiltonian systems would not be possible. Although the gradient performance poor, conventional architectures can still fit functionals and the MSE for validation samples within the training distribution (Base) is low. When considering generalization performance, IKFs can readily generalize to unseen input functions as well as unseen discretizations in the linear case. This empirically observed generalization is theoretically motivated; the Riesz
5

---

A
Linear Functional Derivatives: diFiddu x²
Nonlinear Functional Derivatives: dF /du 3u
1.0
IKF
After training to regress Flu], the gradients of the network with respect to an imput function are plotted for different models. Left, A: Neural functionals are able to accurately model derivatives, even on unseen discretizations, whereas other architectures fail even on their training distribution. Right, B: Nonlinear functionals have more complex derivatives that change with input functions u. In this instance, quartic function is learned and the model extrapolates to unseen discretizations.
Figure 2:
representation theorem guarantees that single kernel Ki can represent a linear functional Fi for all possible input functions u(z) (e.g., polynomial, exponential, etc.). Generalization performance with nonlinear functionals is more challenging, as larger polynomial coefficients c; and larger coordinates I, quickly increase Fu [u[](())] due to large exponents. Despite this, IKFs still outperform other architectures in these domains, and learn accurate nonlinear functionals and their derivatives for samples within distribution.
4.2 1D Advection and KdV Equations
To evaluate the Hamiltonian Neural Solver (HNS) on PDE problems, we consider the ID Advection and the ID Korteweg-De Vries (KdV) equation. The Hamiltonians for each PDE are given as [2]:
Formulation
1 Hasto [u] = 1-000-1---- (u()))
8uu Hadd [u] = - (u(=))) u(x)
(1 dr,
(8)
Ix
Given an initial state u., the Hamiltonian for each PDE is conserved over time (i.e., H(w") = H [u (')). For both PDEs, the operator J is defined as 0g. We can check that Hamilton's 's equations hold for the ID Advection equation: 8Hade uu(x),
J( 8'Had
Ou
des
8u
8u
Ox
which recovers the correct PDE for ID Advection. A similar calculation can be done to recover the PDE for the ID KdV equation (ut + utk + Uezr = 0) from Hamilton's equations by using ou - ugz. Therefore, given a training dataset of spatiotemporal PDE data, the labels can be computed using finite differences at every timestep t. The gradients of the HNS network Va(Ho(x, u²)) are then fitted to these labels.
sti
Evaluation
Data for the 1D Advection and KdV equations generated from numerical solver using random initial conditions 23 24]:
a (x))) = Σ A₁ sin
2xl x
+ 6.)
(10)
L
1=1
where L is the length of the domain, A; E 1-0.5, 0.5). 4 E (1, and di E [0, 2T]: J 5 for Advection and J = 10 for KdV. The Advection equation is solved with n, = 128 and the KdV equation is solved on a domain 2; E [0, 100) with n, 256. Both equations have periodic boundary conditions. For the Advection equation, 1024/256 samples are generated for train/validation, and for the KdV equation, 2048/256 samples are generated. To evaluate temporal extrapolation, the training dataset is truncated to a shorter time horizon. For the Advection equation, it is solved from /. = Os to / 4s during training (nt 200) and evaluated from / = Os to 20s during validation (n; 1000); a similar scenario constructed for the KdV equation with t Os to f 25s for training (n, 50) and t Os to 100s for validation (n, 200).
uniform domain I; [0, 16]

---

M((s over time, Advection
24 (uj over time, KdV
Adv
KdV
110
Metric:
Roll. Err.
Corr. Time
FNO
0.83
68.0
Mis
R
Unet
0.40
134.01106
2
FNO du
0.008
77.61:.6
CN
de
Unet(
at
du
1.057
113.3
Time
501
Time
HNS
0.0039=0.0008
151.1#30
True FNO
Unet
HNS
Table 2: Results for 1D) PDEs. Parame- ter counts are: FNO (65K), Unet (65K), HNS (32K) for Adv, and FNO (135K),
Figure 3: Hamiltonians for Adv and KdV for predicted trajectories, plotted over time. FNO and Unet models shown with their du dt variants. HNS can better conserve Hamiltonians and remain stable over time.
Unet (146K), HNS (87K) for KdV.
Two baselines are considered, a FNO 13 and Unet 25 model. Additionally, the HNS model is used with a nonlinear, FiLM-coodiiionee SIREN kernel. Global conditioning is used in the KdV
experiments by calculating the scale/shift parameters as 7 =[30 (u')]: ]s and 3 [3 (u')];: the global FiLM networks are each parameterized by a shallow, 1D CNN. We find that global conditioning is needed in more complex PDEs since the kernel output = ko (-) can depend on more than just Mi given point z;. In particular, the Hamiltonian for the KdV equation has Uzz term that requires non-local knowledge of u'. Lastly, to ensure fair comparison, we also train variants of the FNO and Unet baselines that predict temporal derivatives du/d: dt (26), as opposed to the usual variants that predict u 1; these use the same numerical integrator during inference as the HNS model.
The metrics that are reported are the rollout enror = ΣT= C(u', ub), and the correlation time and the correlation criterion the Pearson correlation. Rollout error used measure the accuracy of the predicted trajectory; for chaotic dynamics, rollout error is often skewed after autoregressive drift, therefore correlation time measures the portion of the trajectory that is accurate. Each metric is averaged over the validation set and the mean and standard deviation over six seeds is reported in Table Furthermore, given a predicted trajectory [ug..... (uj....3 ug from a model, we can examine its Hami ionian at each timestep and plot the scalar value H[ug] over time, shown in Figure
max(t) s.t. C(u', ut)) 0.8. The loss used is relative L2 loss, or L(u', ut)) -
We observe that the proposed HNS model is able to outperform current baselines, both on quantitative metrics and when conserving Hamiltonians over time. HNS is also efficient, capable of achieving better performance with nearly half the parameters of baseline models due sharing weights between input points Fi, Furthermore, despite training on a shorter time horizon, HNS can extrapolate and remain stable in longer prediction horizons. This more pronounced in the KdV equation, where initial states are smoother than later, more chaotic states. An example trajectory is shown in Figure 5 where the baselines struggle to predict later states since these are not seen during training.
4.3 2D Shallow Water Equations
Formulation To evaluate HNS on a more complex system, we consider the 2D Shallow Water Equations (SWF). These equations are approximation of the Navier-Stokes equations when the horizontal length scale is much larger than the vertical length scale. We consider its conservative form, which is often used simple model for atmosphere ocean dynamics.
For vector-valued functions u(x, y) E Rd, the operator J becomes a d matrix of operators. For the 2D Shallow Water Equations, the operator matrix J is given by (2729
9
0
Our Ou Oy
J
g 0
0
(11)
0
Ör
where q is the vorticity. The Hamiltonian for this system is:
H(u] Hlu Uy M) 1 Mo + -
1 g
(12)
where g is the gravitational constant and the velocities and height 2y Ny h are usually defined on grid (i.c., Che (2,, ), One can verify Hamilton's equations and that J 5Hl
Su recovers the 2D
7

---

Model
Sines
Pulse
NV us) over time, SWE-Sines
7f(me over time. SWE-Pube
Transolver
0.084=000081 0.12240.0074
1.2
FNO
0,057++0002
0.117=0.0009
PINO
0.053++0.005
0.114=0.0003
Mju
1004
Unet
0.010.10.004
0.002.0..006
0.5
HNS
0.026:10.003
0.021 0000015
L,00
0.7
Table 3: Rollout errors for 2D SWE. Parameter counts are: FNO/PINO (7M), Transolver (4M), Unet (3M), HNS (3M). Models are evaluated on in-distribution ICs (Sines) or out-of- distribution ICs (Pulsc).
Time sbop
me stop
FNO
PIND
Tratscheer
Figure 4: Hamiltonians for Sines and Pulse ICs, plotted over time. Despite potentially larger rollout errors, HNS can better conserve the Hamiltonian of predicted solutions.
Shallow Water Equations:
0,h V (vh) 0 Vv= -
(13)
Evaluation Data for the 2D SWE system is generated from PyClaw 30 32 soidal initial conditions (Sines) for the height:
using random sinu-
2mm; Ln
2mm
(x,,)) Σ A sin(
+ Φi)
(14)
Ly
with initial velocities set to zero. Constants are randomly sampled from A; E [-011, 0.1]. my, ny E (1,2, and di [0, 2T] with J H 1, and L, = Ly The domain is (æ y1) E [-1 with ny = 128; each sample has n, = 100 timesteps from t = Os to t 2s. The gravitational constant set to g 1: 256 samples are generated for training and 64 samples are generated for validation. Periodic boundary conditions are used. To evaluate model generalization, we generate an additional test of 64 samples with a random Gaussian pulse as initial conditions (Pulse). The initial height is set to an isotropic Gaussian centered at the origin with random covariance a E [0.1, 0.5] and the height scaled to be between [1, 1.5]; initial velocities are zero. The domain, BCs, and discretizations remain the same.
We consider benchmarking against Unet and FNO models, as well PINO 33 and Transolver 34 models to compare against a physics-infornee method and a newer attention-based surrogate. We instantiate the HNS model with global, FiLM-conditioned SIREN kernel. Validation rollout errors for in in-distribution ICs (Sines) and out-of-distribution ICs (Pulse) are given in Table where means and standard deviations are calculated over three seeds. Additionally, the Hamiltonian over time for model rollouts are plotted in Figure 4 The derivative variants of baselines did not train stably, possibly due to sharp gradients that develop in the absence of viscosity, and are reported. Additionally. the baselines are trained with the benefit of normalization, as the height channel has different scale than the velocity channels; however, since the Hamiltonian not invariant under scaling and shifting, normalization was not performed for the HNS model.
For test samples within the training distribution (Sines), conventional neural surrogates are able to perform well and achieve low loss, and HNS is on-par with these models. Adding a physics- informed Joss in PINO only marginally increases the performance of FNO, which is consistent with its interpretation as a soft regularizer. When looking at the Hamiltonian of predicted solutions, we find that HNS can predict solutions that better conserve physical invariances. We validate this in Appendix by calculating the error of the Hamiltonian of predicted trajectories across different experiments. This ability to conserve the Hamiltonian helps HNS to generalize to unseen initial conditions, and it achieves better performance on the Gaussian Pulse test set. Interestingly. it also conserves the Hamiltonian on out-of-distribution inputs, suggesting a strong inductive bias in the model. To understand what mechanisms enable energy conservation, we also ablate different methods within HNS to understand their contributions, shown in Appendix [.. Lastly, to visualize model performance, predicted trajectories of Sines and Pulse ICs are shown in Figures 0 and

---

Related Works
Hamiltonian Neural Networks Original work on Hamiltonian Neural Networks (HNNs) 35
36
established them as promising architecture to respect physical laws, while later work improved their performance [7 41, investigated the mechanisms behind HNNS [42] or used HNNs in different
A closely related set of works develop neural networks in Lagrangian framework 46 which also preserve learned energies. More broadly, learning an energy and optimizing its gradient is an established framework in ML-based molecular dynamics simulation [48- 52.. Despite research in this area, no prior works have considered learning with infinite-dimensional Hamiltonian systems and most HNNs have only been demonstrated in simple, particle-based systems.
applications 43 45
Functional Approximation While uncommon, there are prior works that investigate learning infiniteddimensional to scalar/vector maps. Most methods are based a kernel integral and are used to functionals (2). construct functional autoencoders 4, model parameler-t0-obseraabb
maps (3. build infiniteddimenionna discriminators for generative modeling [5]. The use of
kernel integral to model a functional is also well studied in machine learning for DFT, since the
Hohenberg- Kohn Theorems guarantee that the total energy is the integral of unknown function of electron density (53)55) Lastly, a separate approach considers modeling functionals through the cylindrical approximation, rather than kernel integral [56].
Neural PDE Solvers Neural PDE solvers are a growing field, with many works proposing ar- chitectures to improve model accuracy [57 34 58) generalizability [19 stability [62 65). Specific to our work, many PDE solvers have also employed physies-based prior knowledge to improve model performance and adherence governing laws. The most straightforward approach is to use physics-informed losses |66 [3 67], but another, more fundamental approach is to enforce physical symmetries and equivariances in the network itself (68/72). Interestingly, by Noether's theorem, every symmetry has an associated conservation law, which HNS relies on as an inductive bias; through this lens the current work is another way of softly enforcing symmetries in predictions.
6 Conclusion
Limitations Despite an interesting theoretical basis and its novelty, Hamiltonian Neural Solvers have limitations that we hope future work can address. Firstly, the theory for nonlinear functionals is underdeveloped. There is also additional overhead during inference to backpropagate through the network and evaluate the linear operator J. Scalability is also a concern due to the reliance on numerical quadratures. This increase in runtime is quantified in Appendix B.2
Additionally, implementing J without neural operator may require care to avoid numerical errors; this is further discussed in Appendix D.I Furthermore, HNS is only applicable to Hamiltonian systems; this would exclude systems with dissipation, although modifications exist extend Hamil- tonian mechanics to con-conservvttee systems 73 776. Even for conservative PDEs, implementing a Hamiltonian structure not as straightforward as training a neural solver with next-step prediction. Lastly, while not limitation, an interesting observation is that introducing loss on H(u') H[us] harms HNS performance; we hypothesize that since the Hamiltonian is conserved over time, learning the identity mapping is the easiest way to minimize this error and thus degrades model training.
Outlook this work, we have proposed novel method for designing neural PDE solvers that respect conservation laws. We verify that kernel integrals are able to implicitly learn functional derivatives, as well as propose parametrization using neural fields. Using the resulting architecture in Hamiltonian framework allows stable and energy-conserving predictions of 1D and 2D PDEs. These capabilities enable Hamiltonian Neural Solvers to generalize to certain out-of-distribution inputs, such as longer time horizon or unseen initial conditions.
While learning functions and operators are dominant paradigms, approximating functionals may also have interesting uses. Functionals provide concise descriptions of physical systems by integrating physical variables, such as describing the energy concise, functional derivatives are often more relevant and can contain dynamical information or perhaps be used for optimization. We hope that future work can continue to improve functional approximators as well as take advantage of these models in new and insightful applications.
molecule or the drag of an airfoil. Although

---

Acknowledgments and Disclosure of Funding
We would like to thank Dr, Amit Acharya for his insightful discussions and help in conceptualizing this work.
References
[1] Jerrold E. Marsden and Tudor S. Ratiu. Introduction to Mechanics and Symmetry. Springer
New York, NY, 1999.
[2] Zheyuan Hu, Tianbo Li, Zekun Shi, Kunhao Zheng, Giovanni Vignale, Kenji Kawaguchi, Shuicheng YAN, and Min Lin. Neural integral functionals. ICLR 2023 Workshop on Physics for Machine Learning, 2023. URL httpss://pppnrreieerttt
[3] Daniel Zhengyu Huang, Nicholas H. Nelsen, and Margaret Trautner. An operator learning perspective on parameter-to-observaae maps. Foundations of Data Science, 7(1):163-225, 2025. ISSN 2639-8001. doi: 10.3934/fods.2024037. URL / //x doi 100/10.0000 3934/ fods 2024037
[4] Justin Bunker, Mark Girolami, Hefin Lambley, Andrew M. Stuart, and T. Sullivan. Autoen- coders in function space, 2025, URL https://arxiv..rgg.. a200 01362
[5] Md Ashiqur Rahman, Manuel A Florez, Anima Anandkumar, Zachary E Ross, and Kamyar Azizzadenesheli. Generative adversarial neural operators. Transactions on Machine Learning Re- search, 2022. ISSN 2835-8856.
[6] E. Kreyszig. Introductory Functional Analysis with Applications. Wiley Classics Library. Wiley, 1991. ISBN 9780471504597.
URL. https://bbokk: http:/boooo..ogglleee../ google..o....oooooo com/bookssid-
AQEMEAAAQBAJ
[7] Allan Pinkus. Approximation theory of the mlp model in neural networks. Acta Numerica, 8: 143-195, 1999. doi: 10.10077509922422000099999
[8] Moshe Leshno, Vladimir Ya. Lin, Allan Pinkus, and Shimon Schocken. Multilayer feed- forward networks with nonpolynomial activation function can approximate any function. Neural Networks, 6(6)8861-867, 1993. ISSN 0893-6080. doi: https://doi.org/10..u
90893-600800058801115. URL https: //wtp: sciencedirect. com ht////icicccarrccllcc p11/$0893608005801315
[9] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Arizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Jourmal of Machine Learning Research, 24(89):1-99) 97, 2023. org/
http: http://mir..rrr/pppprr////////// Dapers//24/21-145
html
[10) Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations, 2020. htt.tttss://arxi http:/ppsss//rrrt.vooo/////////0 2003 03485
[11] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 6755 6766. Curran Associates, Inc., 2020. URL httpps: www.wedddddninnnnnn..eee.eeeee.. wweripp.cc/ppppe.. wwtppcppppppiiieee filee/ pd.ff
[12] Zongyi Li, Nikola Borislavov Kovachki, Chris Choy, Boyi Li, Jean Kossaifi, Shourya Prakash Otta, Mohammad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadeneshelii and Anima Anandkumar. Geometry-informed neural operator for large-scale 3d pdes, 2023.
URL peattp:://axxxvvorr//aaaaa
00583
[13] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations, 2021. URL https: htaps:///axxvv org/ httgg/ab//2000 abs./2000008855 08895
10

---

[14] Yiheng Xie. Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. Neural fields in visual computing and beyond, 2022. URL https://arxiv http:///axxivo..g/bbb..ss/ wttp:bogg////aa..e...oo 2111 11426.
[15] Songyou Peng, Michael Niemeyer, Lars Mescheder, Marc Pollefeys, and Andreas Geiger. Convolutional occupancy networks, 2020. URL htt.:tttss//aaxi.orr//a/a 2ttp//2000004118888
[16] Ishit Mehta, Michaël Gharbi, Connelly Barnes, Eli Shechtman, Ravi Ramamoorthi, and Man- mohan Chandraker. Modulated periodic activations for generalizable local functional represen- tations, 2021. wtt.ltttts::/ariiv..oo org/ abs httg://bbs/200..00m 03960
[17] Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove. Deepsdf: Learning continuous signed distance functions for shape representation, 2019. URL
https://arxiv .rg/ abaaaaas/999 1901. 05103
[18] Lars Mescheder, Michael Occhsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger. Occupancy networks: Learning 3d reconstruction in function space, 2019. URL https
//arxiv.ogg/aat abs/122 03828
[19] Ethan Perez, Florian Strub, Harm de Vries, Vincent Dumoulin, and Aaron Courville. Film: Visual reasoning with general conditioning layer, 2017. URL
1709 07871
[20] Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein. Implicit neural representations with periodic activation functions, 2020. URL
https: attps://arxiv. org/ abs/abb///00 09661
[21] Robert M. Miura. The korteweg devries equation: A survey of results. SIAM Review, 18(3): 412 112-459, 1976. doi: 10.1137/1018076 URL https/////// 10. 1rt..//0//011111/11117777 1137/1018076
[22] Peter D. Miller The korteweg-de vries equation, 2018 URL https::/dtpt......
1sa
W18 Wd8.pdf
[23] Anthony Zhou and Amir Barati Farimani. Masked autoencoders pde learners, 2024. URL
wwtp://r//abb/23.11111 17728
[24] Johannes Brandstetter, Max Welling, and Daniel E. Worrall. point symmetry data augmen- tation for neural pde solvers, 2022. URL
[25] Jayesh K. Gupta and Johannes Brandstetter. Towards multi-spattoooempraa--sclll generalized 15616
pde modeling, 2022.
[26] Anthony Zhou and Amir Barati Farimani. Predicting change, not states: An alternate framework for neural pde surrogates. Computer Methods in Applied Mechanics and Engineering, 441: 117990, June 2025. ISSN 0045-7825. doi: http:/..cmaa225.117799.9. URL http:
doi org/ 1016/j cma. 2025. 117990.
[27] J.S. Allen and Darryl D. Holm. Extended geostrophh hamiltonian models for rotating shallow water motion. Physica D: Nonlinear Phenomena, 98(2):229-22 248. 1996. ISSN 0167-2789. doi: http:://bb.orr:/0.00000001::
http:/111111111111111 011.222888969000002222 URL https: /wwww/wwwww.. sciencedirect. com/ Nonlinear Phenomena in Ocean Dynamics.
[28] Ripa. Conservation laws for primitive equations models with inhomogeneous layers. Geophysical & Astrophysical Fluid Dynamics, 70:85-111, 6 1993. ISSN 10290419, doi: 10. 1L.1880003091992332223333
[29] Paul J. Dellar. Common hamiltonian structure of the shallow water equations with horizontal temperature gradients and magnetic fields. Physics of Fluids, 15:292-297, 2003. ISSN 10706631. doi: 10.10631 10.1063/1.153055
[30] Clawpack Development Team. Clawpack software, 2024. URL http:htt:::///t www.:wwwwwww...... org Version 1.0.

---

[31] David I. Ketcheson, Kyle T. Mandli, Aron J. Ahmadia, Amal Alghamdi, Manuel Quezada de Luna, Matteo Parsani, Matthew Knepley, and Matthew Emmett. PyClaw: Accessible, Exten- sible, Scalable Tools for Wave Propagation Problems. SIAM Journal on Scientific Computing, 34(4):C210-C231, November 2012.
[32] Kyle T Mandli, Aron Ahmadia, Marsha Berger, Donna Calhoun, David L George, Yiannis Hadjimichael, David I Ketcheson, Grady Lemoine, and Randall J Le VeVeqe. Clawpack: building an open source ecosystem for solving hyperbolic pdes. PeerJ Computer Science, 2: e68, 2016. doi: 10.7717/pperj-cs..6
[33] Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informet informed neural operator for learning partial differential equations, 2023. URL URI.ttps:: httts://arrii: org/ ///g//b///// 2111.03799
[34] Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver:
A fast transformer solver for pdes general geometries, 2024. https /tttip:://axxi...o org/ abs/2402. 022//4022.03666.
[35] Sam Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks, 2019. URL httpp::/aaxxv.org/// hat///bb//1900 /ttpp///9....11. 01566
[36] Tom Bertalan, Felix Dietrich, Igor Mezić, and Ioannis G. Kevrekidis. On learning hamiltonian systems from data. Chaos: An Interdisciplinary Journal of Nonlinear Science,
12 2019. ISSN 1054-1500, doi: 10. 10.1063/1.511822) 5128231
(RLLhhtps:://dtt org/ 10. 10.13/101///1111111111
[37] Marco David and Florian Méhats. Symplectic learning for hamiltonian neural networks. Journal of Computational Physics, 494:112495 2023, ISSN 0021-9991. doi: https:/doi..rr/ ic.00///..pp.2223..11 p11/$0021199923065009
URL htttps//wwet scceceeiicccc.....c......... wttp:/reett.ccccciieieeeee.eeeee comm/sciiinc//arricie...
[38] Philipp Horn, Veronica Saz Ulibarrena, Barry Koren, and Simon Portegies Zwart. A generalized framework of neural networks for hamiltonian systems. Journal of Computational Physics, 521:113536, 2025. ISSN 0021-9991. doi: https://doi.or//10000000
2024. 20.2.
URL
https://wwet http://www.ciienccirr...ccc http://wwwaacatc..occccocoocaaaa
[39] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Symplectic ode-net: Learning hamiltonian dynamics with control. In International Conference on Learning Representations,
2020. [RR/httts///pperrviiie
[40] Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. In International Conference on Learning Representations, 2020. URL
https://openreevee: net www.mmrmmmm... Bttp///oooom//dddgPPPPPP
[41] Alvaro Sanchez-Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ode integrators, 2019. URL httttsts///ariiii org abs/a/as///1// 12790
[42] Nate Gruver, Marc Finzi, Samuel Stanton, and Andrew Gordon Wilson. Deconstructing the inductive biases of hamiltonian neural networks, 2022. URL URL_nttps:/arrvvv.or///bbb
04836
[43] Marios Mattheakis, David Sondak, Akshunna S. Dogra, and Pavlos Protopapas. Hamiltonian neural networks for solving equations of motion. Physical Review E, 105(6). June 2022. ISSN 2470-0053, doi: 10. 10.1103/physreve.105.066 105.0 URL //dx: doi. org/10..0 1103/
PhysRevE. 105. 065305
[44] Peter Toth, Danilo J. Rezende, Andrew Jaegle, Séhastien Racanière, Aleksandar Botev, and Irina Higgins. Hamiltonian generative networks. In International Conference on Learning Representations, 2020. U httpttt//ppeeeeeiieeee http:///ppppreeeeeeeeeeet/oororp
[45] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty. Extending lagrangian and hamiltonian neural networks with differentiable contact models, 2021. URI. http:httt//ttt htttps://arxiy org htg/abs/2200..007/
12

---

[46] 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. URL http:p//w/gpppppoorrsssese//e//e
[47] Michael Lutter, Christian Ritter, and Jan Peters. Deep lagrangiann model prior for deep learning, 2019. URL
networks: Using physics as 04490
[48] Stefan Chmiela, Alexandre Tkatchenko, Huziel E. Sauceda, Igor Poltavsky, Kristof T. Schütt, and Klaus-Robert Müller. Machine learning of accurate energy-conserving molecular force fields. Science Advances, 365):e 2017. doi: 10.1 (0..116/sciadd.16600000 URL https //www org doi hwtpoi/aa////i// 1ww. 1101 1603015
[49] Johannes Gasteiger, Janek GroB, and Stephan Günnemann. Directional message passing for molecular graphs, 2022. URL https://arxiv.ortt wwtp//////......./////////////o. 0011/220330011333
[50] Kristof T. Schütt, Pieter-Jan Kindermans, Huziel E. Sauceda, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller. Schnet: A continmuuu-fitterrcoooouuue neural network for modeling quantum interactions, 2017. BRLLattsss//arrvv.og//bb///oo 08566
[51] Linfeng Zhang, Jiequn Han, Han Wang, Wissam A. Saidi, Roberto Car, and Weinan E. End-to- end symmetry preserving inter-atomic potential energy model for finite and extended systems, 2018. TtRL.httts:/// //tp:v///xxvvav/a wtr./ii..//..ooooooooooooooooorr abs/1805.
[52] Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P. Mailoa, Mordechai Kornbluth, Nicola Molinari, Tess E. Smidt, and Boris Kozinsky. E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. Nature Communications, 13(1), Utttttp:://dx.c doi org/
May 2022. ISSN 2041-1723. doi: 10. 10.10888/411477-0222999333 10.11088/411467-022-899 2993999--5
[53] Ryan Pederson, Bhupalee Kalita, and Kieron Burke. Machine learning and density functional theory. Nature Reviews Physics, 4(6):357-358 May 2022. ISSN 2522-5820. doi: 10. 1038/ 112254-022-00477
URt.tttp://dx. doi wtt.:o....//////...//.... 042-2422-2-000000000-022
[54] Tianbo Li, Min Lin, Zheyuan Hu, Kunhao Zheng, Giovanni Vignale, Kenji Kawaguchi, A.H. Castro Neto, Kostya S. Novoselov, and Shuicheng YAN. D4FT: A deep learning approach to kohn-sham density functional theory. In The Eleventh International Conference Learning Representations, 2023. URL http:////p.paeeeeeeeeeeeeeooaaaa
[55] L. Fiedler, K. Shah, M. Bussmann, and A. Cangi. Deep dive into machine learning density functional theory for materials science and chemistry. Phys. Rev. Mater, 6:040301, Apr 2022. doi: 10. 10.1103/hhsRRe Materaal....r............... Physkevaateriills.c 640330
URL htttp:://ind: aps wwwp///////...o/o/////eooooooeet 10.01...... 1103/
[56] Taiki Miyagawa and Takeru Yokota. Physics-informe neural networks for functional dif- ferential equations: Cylindrical approximation and its convergence guarantees. In The Thirty-eighth Annual Conference Neural Information Processing Systems, 2024. URL
https://openrevies wtwp:rrrwwieee.eeeooomuumaaaaeee nett:forumm.dcc552000XC.
[57] Phillip Lippe, Bastiaan S. Veeling, Paris Perdikaris, Richard E. Turner, and Johannes Brand- stetter. Pde-refiner: Achieving accurate long rollouts with neural pde solvers, 2023. URL
https://arxiv.. org/ abs/ 05732
[58] Benedikt Alkin, Andreas Fürst, Simon Schmid, Lukas Gruber, Markus Holzleitner, and Johannes Brandstetter. Universal physics transformers: A framework for efficiently scaling neural operators, 2025. URL http::tttt:::attei org/ abs//422//222//////// 12365
[59] Anthony Zhou, Cooper Lorsung, AmirPouya Hemmasian, and Amir Barati Farimani. Strategies for pretraining neural operators, 2024. URL https://arxiv.. http://axxvvvrx.////// hatg/abs/24000 08473
[60] Zhongkai Hao, Chang Su, Songming Liu, Julius Berner, Chengyang Ying, Hang Su, Anima Anandkumar, Jian Song, and Jun Zhu. Dpot: Auto-regressive denoising operator transformer for large-scale pde pre- pre-training, 2024. URL Ut //arxivv///i///// org/ hts//bbs/4400 abs//403..0. 2454.03542
13

---

[61] Maximilian Herde, Bogdan Raonić, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Em- manuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for pdes, 2024. URL.ttpss: / /arxivv..vvv////////v//aai//iiv org/ /a////b///// ha////bb//20.... 19101
[62] Anthony Zhou, Zijie Li, Michael Schneier, John R Buchanan Jr, and Amir Barati Farimani. Text2pde: Latent diffusion models for accessible physies simulation, 2025, URL https
//axxiv.org www.vv..rg//bbae aa////aa........ 2410.01153
[63] Zongyi Li, Miguel I Huu-Schiaffini, Nikola Kovachki, Burigede Liu, Kamyar Arizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning dissipative dynamics in chaotic systems, 2022. URL htttth::p/aaaatt http://axxxvv..r/bbb.e/
[64] Bjoern List, Li-Wei Chen, Kartik Bali, and Nils Thuerey. Differentiability in unrolled training
of neural physics simulators on transient dynamics. Computer Methods in Applied Mechan- ics and Engineering, 433: 2025. ISSN 0045-7825.
doi: https://doi.or//10..
j.cma 2024. 1144411741111 $0045782524066960
URL hattps:/vee www.:www.eenn...coc www.ssiciiieeeetccooom http:///ci//ccc////rreieee//////
[65] Zijie Li, Anthony Zhou, and Amir Barati Farimani. Generative latent neural pde solver using flow matching, 2025. URL atttss://arae org httg/abb/2503.. 22600
[66] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707,, 2019. ISSN 0021-9991. doi: https://doi.org/100101166 pop. 2018. 10.045. URL
http://www..iieccerrrettt..tt com/
[67] Sifan Wang, Hanwen Wang, and Paris Perdikaris. Learning the solution operator of parametric partial differential equations with physics-informed deeponets, 2021. URL https://arxiv
hrg/abs/2103.1097/
[68] Zakhar Shumaylov, Peter Zaika, James Rowbottom, Ferdia Sherry, Melanie Weber, and Carola- Bibiane Schönlieb. Lie algebra canonicalization: Equivariant neural operators under arbitrary lie groups, 2025. URL https:://tatiiii httpg///aa//21..2//22 02698
[69] Rui Wang, Robin Walters, and Rose Yu. Incorporating symmetry into deep dynamics models for improved generalization. In International Conference on Learning Representations, 2021. URL URRLitttps///operrviiee
[70] Johannes Brandstetter, Rianne van den Berg, Max Welling, and Jayesh K. Gupta. Clifford neural layers for pde modeling, 2023. URL wwt.t:ttts::aiiii wwwp:/ss////a/oee..coeo..osss abs//209.2222.... 04934
[71] Tara Akhound-Sadegh, Laurence Cerreault-Levaseee Johannes Brandstetter, Max Welling, and Siamak Ravanbakhsh. Lie point symmetry and physics informed networks, 2023. URL
https://arxiv. org/ abs 04293
[72] Masanobu Horie and Naoto Mitsume. Physics-embedded neural networks: Graph neural pde solvers with mixed boundary conditions, 2023. URL URL htttps://axxiiii httpp///aarrii.orrr/b/////////
wrg/obs////aa///ooo 2205 11912
[73] Andrew Sosanya and Sam Greydanus. Dissipative hamiltonian neural networks: Learning dissipative and conservative dynamics separately, 2022. URL lattps://axxii> www.ii.or//bb..........o /
2201. 10085]
[74] Harsh Bhatia, Gregory Norgard, Valerio Pascucci, and Peer-Timo Bremer. The helmholtz-hodge decomposition a survey. IEEE Transactions on Visualization and Computer Graphics, 19(8): 1386-1400 2013. doi: 10. 10.1109/VCC22.012....
[75] John W. Sanders, A.C. Deviria, Nathan J. Washuta, Gafar A. Elamin, Kevin L. Skenes, and Joel C. Berlinghieri. A canonical hamiltonian formulation of the ntvier-stokes problem. Journal of Fluid Mechanics, 984:A27, 2024. doi: 10.1017/jfm.200 10.0017/fm..0224.2999
(76] Alexander Figotin and Jeffrey H. Schenker. Hamiltonian structure for dispersive and dis- sipative dynamical systems. Journal of Statistical Physics, 128(4)999--0556, June 2007. ISSN 1572-9613. doi: 10.10007s 1001007777100555500--9911--1 s10955 007799321-1
h /xxx doi. orr.
14

---

[77] Kaushik Bhattacharya, Bamdad Hosseini, Nikola B. Kovachki, and Andrew M. Stuart. Model reduction and neural networks for parametric pdes, 2021. https://///a orw.abssrg//////.
2005 03180
[78] Samuel Lanthaler, Siddhartha Mishra, and George Em Karniadakis. Error estimates for deep- onets: A deep learning framework in infinite dimensions, 2022. (RLttsss:///axx org/
abs/2102. 09618
[79] Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for fourier neural operators. Journal of Machine Learning Research, 22(290): 1-76, 2021.
URL http://jmlr.... wwt.enlrooroooppppepeeeeee wwt..o///pppp/////111.....
html
[80] Taiki Miyagawa and Takeru Yokota. Physics-informed neural networks for functional differ- ential equations: Cylindrical approximation and its convergence guarantees. In A. Globerson, L. Mackey, Belgrave, A. Fan, U. Paquet, Tomczak, and C. Zhang, editors, Advances in
Neural Information Processing Systems, volume 37, pages 72274-7277 72774924409 Curran Associates, neurips cc cc.ppprr.rppppp.ppppppppppp papprrappe_____ files/lpaale pallee/ppaeer/22//fileee 344c77/588rfffff4c00200700006070 Papeer-reeeeee--eeeeeneeeee pdf
Inc., 2024. URL
[81] Daniele Venturi and Alec Dektor. Spectral methods for nonlinear functionals and func- tional differential equations. Research the Mathematical Sciences, 8(2):27, Apr 2021. ISSN 2197-9847. doi: 10.1007/<40067-002--00666 000087-021-0265--11
URL https://doi. 1tt.:oi..rg/11100000..00 10t//10.10011
[82] K. O. Friedrichs and H. N. Shapiro. Integration over hilbert space and outer extensions. Proceedings of the National Academy of Sciences, 43(4):336-338, 1957. doi: 10. 10733/3/pa
4.336. UtLpssts::/ wwwp:///wwwwppp......mmmm pnas. org.//os/abb/l1.1/77//////
43.4.330
[83] Randall J. LeVeque and David L. George. High-resoluiie Finite Volume Methods for the Shallow Water Equations with Bathymetry and Dry States, pages 43-73. doi: 10.1142/ 9789812790910 00002 URL ( http://wwww.rrrlldseeet......... com
htt.//fii.coo/////aaaaaaaaaaaa/a http//doi//a/b///// 1tabb///...11......... 1142/
9789812790910 0002
[84] Nguyen Ba Hoai Linh and Dao Huy Cuong. A well-balanced finite volume scheme based on planar riemann solutions for 2d shallow water equations with bathymetry. Applied Math- ematics and Computation, 457:128167, 2023. ISSN 0096-3003. doi: https://doi.org/I1
12166 URL http:://www..sceeceddirrr. www.coocccccioe wom.//ciinc////tiicc.a...c articcee p11/$009633032300333
[85] Pavel Holoborodko. Smooth noise robust differentiators. 2008.
[86] Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R. Chakravarthy. Uniformly high order accurate essentially non-oscillatory schemes, iii. Journal of Computational Physics, 131(1):3-47, 1997. ISSN 0021 doi: https://doi.
19tp://cccp.199.6666 5632. URL
http://dccectcccoocccc/////////e https://www/ wtt.://ccc..cc........ca...aaaaa
[87] Xu-Dong Liu, Stanley Osher, and Tony Chan. Weighted essentially non-oscillatory schemes. Journal of Computational Physics, 115(1):200-112, 1994. ISSN 0021-9991. doi: https///oi..rr//
10. 10.1006/jcch.119441 1994. D11/500219999884118899
[88] Kumar Rahul and S.N. Bhattacharyya. One-sided finite-difference approximations suit- able for use with richardson extrapolation. Journal of Computational Physics, 219(1): 13-20, 2006. ISSN 0021-9991. doi:
URL
15

---

A
Additional Visualizations
True
FNO
Unet
125
200
200
200
203
140
150
150
150
150
0.75
100
100
0.5
9.00
0.25
0
123 0
Figure 5: Predicted trajectories of the KdV equation. The white dashed line denotes the training 50 are seen during training. Despite not seeing later,
horizon, where only states from more chaotic states, HNS still extrapolates beyond the training horizon by conserving the Hamiltonian, while other models struggle. FNO and Unet models are shown with their du /tt variants.
SWE-Sines, t=10
True
Unet
HWS
125
128
125
125
1.5
64
64
64
0.5
0:00
0
120
128
4
120
SWE Sines, t=50
3.15
1.05
0
120
120
2.50
Figure 6: Predicted trajectories of the Shallow Water Equations (SWE) with sinusoidal ICs, displayed
10 and 50. All models are able achieve good performance on m-distribution samples.
SWE- Pulse, t=10
True
FNO
Unet
HNS
120
329
120
328
1.2
1.50
64
64
1.55
1.52
0
64
120
156
120
136
SWE-Pulse.
L=50
True
FNO
Unet
HINS
1.00
1.00
Figure 7: Predicted SWE trajectories with Gaussian pulse ICs. displayed 10 and 50. Despite showing different behavior, HNS can better generalize to OOD samples.
16

---

B Additional Results
B.1 Hamiltonian Errors and Inductive Bias
Model
Adv
KdV
SWE-Sines
SWE-Pulse
FNO
2.18
NaN
0.0091
0.1326
Unet
0.22
2.37
0.0053
0.0158
FNO( du
2 0.033
6.75e8
Unet( 0.043
5.76e6
HNS
0.0002
1.32
0.0015
0.0003
Table 4: Relative L.2 Error for each experiment, evaluated on the Hamiltonian of predicted trajectories. Despite not including the Hamiltonian in the training loss and testing OOD samples, HNS is exceptional at predicting solutions that conserve the Hamiltonian.
In the main body, we presented qualitative evidence that HNS is able better conserve the Hamil- tonian of predicted trajectories. To quantify this claim, we calculate the Relative L2 Error of the Hamiltonian of predicted trajectories across all PDE experiments. Specifically, the error of sample is given by: # E
H'ulll? where u' is the true solution at time f. up is the predicted Mllill
solution at time I, and HH} is the Hamiltonian given by the integrals in Equations 8 and 12 The Hamiltonian errors are averaged across the validation sets and reported in Table 4 We note that in the KdV experiments, the dynamics are chaotic and the Hamiltonian is highly nonlinear, both of which contribute to large errors once predicted trajectories diverge. Overall, we find that HNS has exceptional capabilities in implicitly conserving energy-like quantities in their prediction.
Interestingly, this capability exists even when not trained with a loss on the Hamiltonian (i.e., ||H[[') 1l33 and when predicting out-t--distribution trajectories. Where does this inductive bias to conserve the Hamiltonian come from? Based on Gruver et al. [2), we hypothesize four inductive biases in HNS that may contribute to this:
1. ODE Bias: HNS predicts du and uses an ODE integrator evolve PDE dynamics. de
2. Hamiltonian Bias: HNS relies on J&H calculate du dr
3. Gradient Learning Bias: HNS relies on autograd(Hslui) u) to calculate 6H Au
4. Neural Functional Bias: HNS relies on integral kernel functionals to calculate (Hty u).
We note that these biases listed from simplest to most complex and are cumulative; for example, is not possible to implement Hamiltonian structure without also implementing an ODE integrator. We modify Unet to incorporate each inductive bias and evaluate its performance on the KdV equation, which is the most chaotic setting. The ODE bias was considered in the main work by training du dt model variants, while incorporating Hamiltonian bias done by changing the training label to 87//u and using J during inference. To extend Unets to gradient learning, we add linear head to project the output field scalar and use automatic differentiation obtain 8H. After training, we evaluate model correlation time and Hamiltonian error on the validation set and the results are reported in 5
We find that the main contributors to HNS performance and its ability to conserve energies over time using Hamiltonian structure as well the neural functional architecture. Using an ODE integrator alone degrades model performance and implementing gradient learning with Unets also harms performance. Interestingly, using gradient learning with conventional architectures disproportionnatly harms energy conservation; this implies that in PDEs, gradient domain learning improves energy conservation only when architectures have meaningful functional derivatives. This is not readily apparent, observations from molecular dynamics or HNNs suggest that gradient domain learning generally improves energy conservation without considering the underlying architectures (35 46 48
Based on these findings, we hypothesize the following mechanism underlying HNS: integral kernel functionals enable the learning of functional derivatives, combining this with difiite-dimeniionall Hamiltonian mechanics leads to its performance and conservation properties.
17

---

Metric:
Correlation Time (t)
Hamiltonian Error (L)
Unet (Base)
125.25
2.37
Unet (ODE)
120.5
5,76e6
Unet (Ham.)
143.75
2.06
Unet (Grad.)
141
4.71
HNS
151.75
1.32
Table 5: Correlation time and Hamiltonian errors for Unet models with increasingly more inductive biases, compared HNS the KdV equation. Using ODE integrators or gradient domain learning both degrade performance, while using Hamiltonian structure or neural functionals both increase performance and energy conservation.
B.2 Timing Experiments
Model (#Params)
Adv
Kdv
SWE-Sines
FNO (65K/135K/MM)
0.083
0.091
0.967
Unet (65K/1466/3M)
0.138
0.146
1.345
HNS (32K/87K/3M)
0.126
0.228
4.547
Numerical Solver
0.000
110.9
33.26
Table 6: Computational cost inference step for each model across different PDEs, given in milliseconds (ms). Model sizes are also reported for the (Adw/dv//WEE) experiments. Each run uses single NVIDIA RTX 6000 Ada GPU. Numerical solvers on AMD Ryzen Threadripper PRO 5975WX 32-Core CPU due to lack of GPU compatibility.
The additional overhead of using automatic differentiation, evaluating the operator J, and applying an ODE integrator during inference increases the computational cost of HNS. We examine this in Table by performing a model rollout and averaging the time per inference step across a batch. Additionally, we report the runtimes of the numerical solvers used. In the case of 1D Advection, the analytical solution used. Randomly sampled initial conditions can cause variability in the stiffness of the solution, therefore, solver runtimes are averaged across 10 samples. We note that the 2D SWE solver (PyClaw) [31 uses a Fortran compiler and is optimized for hyperbolic PDEs, while the KdV solver is purely Python-based. We note that for fairer comparison, the numerical solvers should be coarsened to match the accuracy the neural solvers.
Consistent with prior works, FNO models remain a fast neural PDE solver. While HNS may have lower parameter count, this may differentiation calculates the gradient with respect to each input point, therefore the computational cost of this step can scale with the size of the input grid. This can be seen when comparing the Adv and KdV experiments, where the grid size doubles from nz 128 to = 256, and the computational cost roughly doubles as well. In 2D experiments, the cost of HNS is appreciably more than other baselines, although it is still modestly faster than optimized numerical solver. We do not anticipate that the finite-difference (FD) schemes used in J or ODE integrators contribute meaningfully to computational cost, since FD schemes use (MM) operations, where M is the number of grid points, and numerical integration occurs in constant time.
misleading measure of inference speed. Automatic
B.3 Ablation Studies
There many choices for parameterizing a kernel integral, such as from different kernels (linear, nonlinear), different architectures (MLP, SIREN), different conditioning mechanisms (concat/FiLMM and different receptive fields (local, global). Since this design space is large, consider the effects of incrementally improving each aspect and report the performance different HNS models Table 7 We find that nonlinear kernel is needed in all cases, since all the Hamiltonians tested nonlinear. In addition, the use of FiLM conditioning and SIRENs tend to improve performance. Lastly, global 0
conditioning is necessary for the KdV equation, which has nonlocal terms
in the Hamiltonian. 0g
18

---

Model
Rollout Error
Model
Correlation Time
Model
Rollout Error
Base
1.248
Base
51.75
Base
0.113
Nonlinear
0.028
Nonlinear
73.00
+ Nonllnear
0,034
+ FiLM
0.016
+ FFLM
65.25
+ -FLLM
0,030
-SIREN
0.003
+SIREN
73.25
+SIREN
0.026
+ GGlbbal
0.057
+ Global
151.75
+ GGlobal
0.024
(a) Ablations for ID Advection
(b) Ablations for ID KdV
(c) Ablations for 2D SWE-Sines
Table 7: Ablation studies for the HNS architecture. The base architecture is specified by (Linear, Concat, SIREN, Local). Each successive row is cumulative, adding an additional feature to the model. Parameter sizes are approximately constant across ablations.
B.4 Varying Dataset Sizes
We study the effect of varying the number of samples the training dataset on HNS performance in Table 8) The training dataset is truncated by a factor of 8, 4, or 2, and evaluated on the same validation dataset. This repeated across all experiments/PDDSs In general, increasing the dataset size consistently improves model performance. This is more acute for 1D) equations. For 2D SWE, the extra spatial dimension may allow for more data per sample relative to the complexity of the dynamics.
Dataset Size
Adv
KdV
SWE-Sines
SWE-Pulse
f=1/8
1.320
71.00
0.0293
0.0243
f=1/4
0.302
71.25
0,0274
0.0230
f=1/2
0.075
95.25
0.0254
0.0227
f=1
0.003
151.75
0.0249
0.0216
Table 8: Rollout error (4) or correlation times (f) HNS models trained on datasets downsampled
by a factor of (e.g.. f = 1/2 uses half the number of training samples in the training set).
B.5 Discretization Invariance
We experiment with querying models at different discretizations than its training resolution. The results are Table g with training resolution in bold. Rollout errors are averaged over the validation set, where validation labels are either downsampled or re-solved at higher resolution using the same initial conditions. Additionally, for 1D Advection, we conducted an experiment in which the model is queried
unseen grid z [0, 32], n. 256, after being trained grid I [0, 16], n2 = 128, denoted as 256*.
(a) Rollout Error Different Resolu-(b) Rollout Error at Different Resolu-(c) Rollout Error at Different Resolu- tions, 2D SWE (Sines)
tions, ID Adv
tions, 2D SWE (Pulse)
Resolution
HNF
FNO
Unet
Resolution HNF FNO Unet
Resolution HNF FNO Unet
64
0.0029 0.064 1.19
(64, 64)
0.028 0.066 0.098
(64, 64)
0.035 0.173 0.173
128
0.0033 0.064 0.075
(128, 128)
0.029 0.066 0.008
(128, 128)
0.024 0.174 0.117
256
0.0055 0.050 1.32
(256, 256)
0.026 0.066 0.1
(256, 256)
0.019 0.174 0.225
256*
0.0021 1.34 0.072
Table Rollout error of models at different resolutions for 1D Ady and 2D SWE datasets.
In all cases, HNF's are capable zero-shot super-resolution, as is possible with Neural Operators. Additionally, the error is approximately constant across discretizations and on new, unseen grid for ID Advection. Predictably, FNO models have nearly constant error across discretizations. Interestingly, for extrapolated grid in ID Advection, FNO struggles to make predictions. For Unet
19

---

models, the eπror increases when the grid spacing changes, however for a constant spacing, it is able extrapolate to unseen grids (i.e. [0, 16] [0, 32) with dr unchanged).
C
Additional Theory
C.1 Proof of Linear Functional Approximation
Interestingly, while development of neural functionals was arrived upon independently, a proof of linear functional approximation exists. In fact, approximating linear functional with an integral kernel an intermediate step in operator approximation theorems [77] and its proof comes as a series of lemmas in Kovachki et al. [9]. For the sake of completeness, provide relevant definitions and cite the necessary lemmas here.
Definitions (Appendix A, Kovachki et al. [9] Let C Rd be a domain. Let No = U {0} denote the natural numbers that include zero. Let X be Banach space, and X be its continuous dual space. Specifically, the dual space contains all continuous, linear functionals f : - R with the norm: || SUPpEx,]|æ|| [(xx) < co.
For any multi-index a E Ng. 0 is the a-th partial derivative f. The following spaces are defined for m E No:
C(D) = ( D + R: f is continuous]
Cm((D) [f: : - + 80 E (D) |a|1 s m}
(D) [f (D)) supp(f) c D is compact)
(D) [ (D) maXuglal Sm supzeD f(x) oo}
. CD) U E (T) 0 is uniformly continuous A0 |//| I m}
Additionally, we cite definitions of neural networks from Section 8.1, Kovachki et al. 9. For any and R. neural network with layers is given by:
Na Rd. Rd) (f Rd RRd ff) (. W.o....aoWWoz bo) b....)
+ bn. Wn E Rdd xdn-1 Wo E Rdoxd Wi E Rd xdo
bo 12 b1 Rd ,...,bb, by Rdd do, di...., d' N
The activation functions a are restricted to the set:
a E Am == fo C(R): In E N s.t. N, is dense in CTT(K) C Rd compact]
to allow universal approximation. Lastly, functions restricted to be real-valued although exten- sions to vector-valued functions are possible (Section 8.3, Kovachki et al. (.
Theorem C.1 (Lemma 30, Kovachki et al. 19). Let D C Id be a domain and L E (C(D)) be a linear functional. For any compact set K C C(D) and > there exists a function C such that:
-
S ruudr <e
sup [L(u) aeK
This result establishes that linear functional acting functions in C(D) can be approximated by a smooth integral kernel with arbitrary error e. To show functional approximation with a neural network, slightly modify existing lemma:
Theorem C.2 (Lemma 32, Kovachki et al. (.). Let D C. Rd be a domain and A C(D). Let
E A be a linear functional. For any compact set K C A. a E Ag, and e > 0, there exists number L E N and neural network Ko NN(o; Rdd R) such that:
sup [T(u) ueK
Ko (o(x)u(z)dr)) <c
20

---

Proof. Since K is bounded, there exists number M > 0 such that:
sup wEK
| M
By Lemma 30, there exists function E C (D) such that:
-J
sup uEK
[(u))
kuddr
Since Ag. there exists some L N and neural network E (o; R) such that:
|| xllc S 2MD)
Then for any E K:
[L(u)= - (o(z)u(z)dzh S [(u) - kudæl1 + (N ko)udx|1
18 + M DI * || Kallc
IN
Therefore, linear functionals for continuously differentiable function classes can be approximated to arbitrary accuracy by neural network using an integral kernel. Extensions to Sobolev and Hilbert spaces can also be proved in similar manner.
Interestingly, this overlap between neural functionals and neural operators suggests a deeper connec- tion between the two. The integral kernel functional bears many similarities to the integral kernel operator and, in practice, implemented as an integral kernel operator evaluated at a single output. From an implementation perspective, the two architectures are equivalent in this manner, examining the theoretical basis also suggests that neural functionals are subset of neural operators when approximating linear functionals. This can be seen by looking at operator approximation theorems (Lemma 22, Kovachki al. (9). We refrain from rewriting the theorems here, but focus on necessary modifications. Specifically, to use neural operator to approximate linear functional, the output Banach space y is set to R and intermediate dimensions J, are set 1. The linear map F, X + R from the input field X' becomes F5 wi (z). where wi (z) E X Lastly, the fiite-dimensional map is defined as &p E C(R, IR) and the output map is defined as G,y K +
C.2
Connection to the Cylindrical Approximation
Miyagawa and Yokota (56] propose an alternate mechanism to approximate functionals based on the
cylindrical approximation. We explore possible connection with the current method of using integral kernel.
We consider functionals Flu]: H + I, where we restrict functions u(x) H to be in Hilbert space H with inner product (., 1 ær. Since functions can be represented as a sum of basis functions, u(x) can be written u(z) Σe k=1 (u, (x))x(z). (z). with an orthonormal basis {Φ1 .). Truncating the sum to finite number of terms m allows the cylindrical approximation functionals [80
This is obtained by substituting the truncated series into the functional:
f(a1, a2,...., am)
F
(+))(=))]
(15)
where are basis coefficients (u, Φx). Using the cylindrical approximation f. functional derivatives can be approximated by:
8F
m
0f
2 (x)
(16)
Su
Day
k=1
21

---

which converges to the true functional derivative as 0c 81I When discretizing the functional
derivative, we define a set of n points x [2 ,g,..., Fn at which the basis functions are evaluated, which results in the discrete approximation of the functional derivative:
of ox Bap Σ Bay
m 0f Σ A=1
m Of 3ax
on (222)
(z)
(x(x))
(17)
Bax
Σ
85 Φu Bux
We are interested if the gradient of neural functional that is parameterized by a kernel integral can represent this discrete object, and what conditions needed do Furthermore, there are two levels of approximation (a truncated basis and discretized domain) and we are interested if the gradient of a neural functional can recover functional derivative in its limit. We will work through three cases of increasing complexity, using a linear kernel, a nonlinear kernel, and a nonlinear and global kernel.
Linear Kernel Integrals
A functional parameterized by a linear kernel integral can be written as:
Folu] -
(oooooz)uee (o(z))(z)ddz Ko ko(x;)u;
(18)
using a Riemann approximation and quadrature weights μ4. Without loss of generality, assume some constant quadrature weight μ and write the gradient of the Riemann sum as:
8hg (No(x1))u) (me (x2(())))) (22))))
(x)(11)) ) No (x2))
Sug
VV
Ko ko(z;)u;)=
p
(19)
a Jun
(ma (I )u,)))
1so (xn))
We can interpret this gradient as the discretization of a single basis function Kg Φ1, and without the necessary coefficients af Evidently, this is not enough expressivity approach day"
of m
coefficients and basis functions, even as + co. Intuitively, the output F& [u] does depend on enough, causing its gradient V to be too simple. A potential remedy is to modify the kernel to be nonlinear by using (a, u())).
Nonlinear Kernel Integrals
Following the same approximation, the gradient of a nonlinear kernel
integral can be written as:
a (no u1 u1) dug (no(x2, U2 u2)
2 Kg (Z1, u1) + 15 (0 u1) W1 Byg
Syn
U2 Buz Ko (2 U2) + 10 u2)
VV
Ko (o;; unuus) =
0 Bu.
(I un un))
9 Un Dua
Ky (Ins un) + (xn un).
(20)
using the product rule. In general, Ko Ka(x;; un) be quite complex and nonlinear, however, by construction it only depends on (æi, uz). We define function dxxx)) Du; 1o0 (22,, uz) lu which evaluates the partial derivative 24. This lets us construct most n basis functions {d1, (di,d2,... .g...., d, } however, this is still not enough expressivity, as each entry of the discretized functional derivative still restricted to one basis function, albeit a different one at each row. Despite this shortcoming, writing this expansion allows us to see that producing sum over multiple basis functions requires the kernel to depend globally on u [ui, #2, .2...., Un], since the gradient of the Riemann sum no longer degenerates to single term at each row.
22

---

Nonlinear Global Kernel Integrals Using the same approximation, the gradient of a nonlinear, global kernel integral can written as:
Σ Ke (2 u)us)
Tuu
Ki (2ig u)us)
V.
Ko (2 u)u;) =
Bug
(21)
Kew (zs, u)u;)
The expression ouz
a
Kg ( u)w:) can be expanded into:
a
a
Kig (o u)u;) = 41
ni (11, u) +...
(22)
Jug
+ uj Ou,
ke (2j, u) + ke( u) +...
+ Un Ouj
k0( u)
In general, we can construct a kernel Ki such that simple example is provided based on FiLM modulation:
Bu, Ko u) is different for i = (I.....n.. A
(T( u) = [W" (W x; + + + [b"};
(23)
where W E Rax", W E g1 6" E I2"" 6" I² and [-]a is an operation that takes the i-th row of a column vector. this case, Bu, Ki (x; u) is equivalent to function dzy(x;) W W* With even deeper and more complex kernels, functions dij can become more expressive. Without loss of generality, can also absorb the basis coefficients into dej.
With this formulation, we can see that the j-th row of the gradient vector can be written as a sum of n functions Ei=dy ddy(2;), which reduces to the cylindrical approximation when letting dıy be arbitrary. Furthermore, letting n + oo increases the number of functions dij and coordinates I;, and consistent with the interpretation of simultaneously increasing the number of bases to approach the true functional derivative as well as converting from a discrete to continuous representation.
The cylindrical approximation makes assumption the linearity of F, therefore, when using a global, nonlinear kernel, kernel integral can represent a functional derivative. This theoretical insight is also supported by empirical observations. Through ablation studies in Section [B3] we find that modeling complex functional derivatives (such as the KdV equation) only possible with global, nonlinear kernel, and simplifications of the kernel degrade the performance.
C.3 Derivations with Varying Coefficients and Source Terms
Using HNS with PDEs with varying coefficients and source terms is possible on a case-by-case basis, requiring additional effort to derive the Hamiltonian structure. We provide examples below by introducing coefficients to ID Advection and ID KdV as well as varying bathymetry (source terms) to 2D SWE.
Varying Coefficients In the presence of coefficients, the Hamiltonians for ID Advection and ID KdV become:
Hadv [u] (((=))
Hedulu [u] =/6 (u(x))² (u)))))))) -u((x) (z)dr,
8u (24)
Oz
For ID Advection, one can verify the Hamiltonian structure:
8Hadv
ccu(x), J 88Had -
Ou du Ix di
(25)
8u
du
similar calculation recovers the PDE for the 1D KdV equation (u, + oanz + Buzzz = 0) from Hamilton's 's equations by using &Had &u
w² Bu,a. accommodate varying coefficients, this suggests training HNS with cocfficient information using the loss: L(V He c), 6
du ). In special
cases where the coefficients can factored out (i.e., 6H[m.d] 6H([u ), once can scale the base du
model output Ho x) by to achieve the same effect.
23

---

Varying Source Terms The 2D shallow- watter equations describe free surfaces where the horizontal length scale is significantly larger than the vertical length scale. Within this setup, the free surface can be above a spatially-varyy the elevation of the bathymetry can accelerate or decelerate flows, which acts as source term on the momentum balance
bathymetry b(x; v), or the topography of the bottom surface. Changes in
88 84).
In particular, the shallow-war waater equations are modified to become:
0gh + V (vh) - O,
vvv Vv = -gVh - ghVb
(26)
Furthermore, h(z, y, now represents the height above the bathymetry b(z, y, t). This source term is represented in the Hamiltonian as:
H [u] [M By hl 1 102 n(u M0 +
1 (1+
(27)
and the operator matrix J is left unchanged (Equation TT]. equations recovers the modified shallow-water equations. When comparing this derivation one with constant bathymetry, one can show that:
One can check that applying Hamilton's 's
8Ml [u] 8u
8Huu 8u
0 0 Light
(28)
This suggests that instead of scaling model outputs V. Ho (u, x) (as in the case of coefficients), source terms can be potentially captured by adding terms to model outputs. In PDEs where a decomposition like this is not possible, the source term will need to be added to the model input produce the output VoHoou, x, b), where bold terms indicate discretized fields/coordinates.
D Implementation Details
D.1 Numerical Methods
Numerical Integration ODE integration is performed using a 2nd-order Maass-Basffotth scheme. The solution the next timestep u + is calculated using the current timestep u' and an estimate of the current derivative du di lt= f(u'):
u u² + Al 3 (2 )
(29)
This can be implemented in constant time by caching prior derivative estimates {f(u') : i t). allowing the ODE integrator to have 2nd-order truncation error (O(At)) while remaining fast. At the first timestep. there are no cached derivatives, therefore first-order Euler integrator is used: u = u' + (u'(u'). For the given experiments, these integrators are accurate since At is small enough; to be used with larger timesteps At, higher-order integrators may be needed or smaller, intermediate steps will need to be taken (26). This adds additional consideration to the training and deployment of derivative-bast neural PDE solvers.
Hamiltonian structure involves approximating the operator J. either through numerical methods or neural operators. We opt for numerical methods, as the operators are simple and in all our tests only involve first-order spatial derivatives o,,, For the ID Advection equation, we implement a 2nd-order central difference scheme. For the 1D KdV equation, the spatial derivatives exhibit larger gradients, therefore we use 8th-order central difference scheme with a smoothing stencil to approximate O, [85]. A similar scheme is implemented in 2D approximate 0 and a, for the 2D SWE experiments, and additionally, Savitzky-Gola Golay filter is used to further damp numerical oscillations. These arise since the lack of viscosity creates discontinuities in the shallow-wad flux-limiting or non coniiiillattory schemes [86] 87]. While these implementations work, they still required tuning, as such, the use of neural operator approximate to make HNS fully learnable would be an interesting future direction.
Finite Difference Schemes Implementing
equations, alternative methods may include
Additionally, to calculate temporal derivatives du/di during training, 4-th order Richardson ex- trapolation is used to provide accurate derivatives. the boundaries 0, T a one-sided Richardson extrapolation is used (88).
24

---

Numerical Solvers The analytical solution to the Advection equation is used to generate its datasets. The numerical solver for the KdV equation is from Brandstetter al. [24]: due to the stiffness of the PDE, the solutions do not always conserve the Hamiltonian, 'hee WENO schemes used allow for more stable derivatives, but the cost of numerical viscosity, which can affect the Hamiltonian since it is highly nonlinear. Therefore, we filter the generated KdV solutions and extract 50% of the trajectories with the smallest variation in the Hamiltonian. Similar fluctuations in the Hamiltonian are seen in SWE trajectories generated from PyClaw (31, but the variations are small enough to be ignored.
D.2 Experimental Details
Hyperparameters Each model has different hyperparameters for different PDEs (Adv/KdV/SWE). We provide key hyperparameters for these experiments in Tables TO
12
Adv
KdV
SWE
Adv
KdV
SWE
Adv
KdV
SWE
Modes
10
12
12
Width
8
12
32
Width 32
32
64
Width
24
48
Bottleneck Dim
32
48
128
Kernel local
global global
Layers
5
5
5
Layers
6
6
Cond. FiLM FiLM FiLM
Table 10: FNO Parameters
Table 11: Unet Parameters
Table 12: HNS Parameters
Computational Resources All experiments were run on a single NVIDIA RTX 6000 Ada GPU; ID experiments were usually run within 10 minutes, while 2D experiments were usually run within an hour. The largest dataset around 20 GB, for the 2D Shallow Water equations.
25

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The authors were careful not claim much in the abstract and introduction, focusing on summarizing the paper and providing an introduction to the field and how the paper fits in.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
Limitations are discussed in Section ?? (Discussion).
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
26

---

Answer: [Yes]
Justification: An informal proof is provided in the main body in Section 3, but a formal proof is given in Appendix C.I
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proof' in the paper should be numbered and cross referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have tried to provide the necessary detail for reproducing results, and all code and datasets will be also be released.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reprodueibility may be tricky in some cases, in which authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
27

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
All code and datasets will be made publicly available.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips //ips.cc//pill wttp:ccppppcccceeeeeeeeee..c.ett httpa//gggeeee/eeeuuebiisssssss /Cc////oddiissssseeeeeee ssiiiiiiccyyyyyyy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification: We have provided information on the experimental setup and specific details can be found in the config files of the provided code.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Error bars are given for the main quantitative, PDE results in Tables and Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi-
dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
28

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
Details on the compute resources are provided in the Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification:
The paper does not violate the code of ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- cration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The work is largely foundational and in the field of physics; at this moment there are no likely negative impacts of the work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
29

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
There are no risks for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The authors have created all assets used in the paper.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in package should be provided. For popular datasets, paperswithcone com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
30

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer:
[Yes]
Justification: New code and datasets will have accompanying documentation. Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The research does not involve crowdsourcing or human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Bthies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing or human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, not include any information that would break anonymity (if applicable), such the institution conducting the review.
31

---

16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: LLMs were not used in the core method development of the research.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to L.LM policy for what should or should not be described.
32


Reviews:
[{"review_id": "HeEOPduRIF", "reviewer": "Reviewer_qFJp", "summary": "This paper introduces a neural network ansatz for learning the Hamiltonian functional, which can be used to model dynamics governed by a PDE. The paper parameterizes the functional as a local integral functional, whose functional derivatives can be obtained via AD. The approach is evaluated on 1D advection, KdV, and 2D shallow water equations, demonstrating improved energy conservation and stability compared to existing neural PDE surrogates.", "questions": "- In line 101, you mention that there exists $f(u)$ that can serve as a Riesz-like representation for a nonlinear functional. However, I don't see that there is any guarantee that $f$ can be approximated with $\\kappa_{\\theta}(x, u(x))$. For example, if $H[u]$ is non-local, then $\\kappa_{\\theta}(x, u(x))$ can never approximate $f$ since it only depends on $u$ pointwise.\n- Why is it that in algorithm 1, you train on the functional derivative, while in the toy example, you train on the kernel? It is confusing that you also mentioned that fitting to the kernel also yields the correct functional derivative for the toy problem. If that's true, why change it for the actual problems?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n- The mathematical background on Hamiltonian mechanics is appropriately introduced so that the paper is self-contained\n- The author tested the OOD performance of the model\n- Compared to previous methods, the proposed method captures the PDE dynamics in a more principled way\n- Demonstrates practical benefits, including improved stability and energy conservation, especially for long rollouts\n\nWeaknesses\nThere are existing works (e.g., https://openreview.net/pdf?id=aQuqw6eVKP) that have essentially the same setup, albeit even more general (semi-local functional, whereas in this paper the functional is local). The novelty of the paper comes from applying neural functionals to the specific setting of Hamiltonian mechanics. I feel that the author should clearly state the paper's contribution both in the title and the abstract. The keyword Hamiltonian should be mentioned in the title. Personally, I had very different expectations before reading the paper, by looking on the title.\n- Using quadrature to approximate the integral makes this method not scalable in the domain dimension\n- The reason that the functional derivative here can be directly be computed directly with AD could be explained in a much clearer way, for example, using the Euler-Lagrange formula\n- Typo on line 144", "was_revised": true, "final_justification": "I will maintain my initial score, since the author has resolved my concerns and questions during rebuttal, and the rebuttal with other reviewers did not bring up new information that would overrule my initial assessment.", "rebuttal": "Thank you for taking the time to review our work and provide feedback. Please find a response to your questions/concerns below:\n### Related Works and Framing\nThank you for providing the existing work. Since the submission, we have also found other works that propose functionals approximated by an integral kernel [1,2], and other reviewers have also pointed out similar works [3, 4]. We arrived upon the idea independently, which unfortunately accounts for much of the language/framing used in the title and paper. We are happy to change this, especially since we have a novel perspective on using neural functionals in Hamiltonian mechanics to evolve PDEs. Indeed, this motivates much of our paper, such as verifying the validity of functional derivatives on toy problems and investigating better architectures (SIREN, FiLM). The title will be changed to: “Hamiltonian Neural Solvers through Functional Approximation” (or something similar), and the paper will be revised accordingly. In particular, the related works will be expanded to include prior work on integral kernel functionals and Section 3 will present intuition on their capabilities from a theoretical standpoint, rather than a derivation. Any language suggesting to be the first work on approximating functionals will be changed.\n### Quadratures and Scalability\nThis is a limitation, especially in 3D where numerical integration is very expensive. One potential avenue to alleviate this is to train on low-resolution data and only use the desired resolution during inference. We find that HNFs can be discretization invariant, which we verify in experiments suggested by Reviewer fiA5. \n### Euler Lagrange formula\nThank you for bringing this to our attention, it is a very interesting perspective. We will discuss the relationship of the Euler Lagrange equation to functional derivatives, as well as go through some examples, such as with the 1D Advection equation. \n### Approximation Guarantees\nYou are correct that an assumption of compact support would be necessary to approximate a functional with an integral kernel. In PDE problems this is usually a valid assumption, as most practitioners are interested in a physical response in a predetermined domain. Beyond compact support, we will make it clear that more assumptions are needed about H[u], as it must exist in the dual space of functions that are: (1) continuous, (2) have a uniformly continuous derivative, (3) have compact support, (4) have bounded derivatives. These are given in Appendix C.1. Similar restrictions also exist for Neural Operators in order to achieve their approximation guarantees. \n### Training Differences\nThis is a great question, and one we will clarify in the work. When constructing the toy problems, we tried to use Algorithm 1 to train the FNO/MLP models, however they did not succeed presumably since their gradients calculated by AD don’t allow them to learn underlying functionals. Therefore, to make things fair in the toy problems, we trained all models to fit the scalar-valued functionals, which has the added benefit of evaluating if models can implicitly learn functional derivatives. For the rest of the work, we also tried fitting the scalar-valued Hamiltonians, however, this had an unintended consequence. Since the Hamiltonian is conserved over time, models would usually learn the identity map for its functional derivative, resulting in the PDE not evolving in time (but still conserving the Hamiltonian). To allow the model to learn dynamical information, fitting it to the functional derivative was necessary. \n\n### References\n1. J. Bunker et. al. Autoencoders in Function Space, https://arxiv.org/abs/2408.01362\n2. D. Huang et. al. An operator learning perspective on parameter-to-observable maps, https://arxiv.org/abs/2402.06031\n3. M. Rahman et. al. Generative Adversarial Neural Operators, https://arxiv.org/abs/2205.03017\n4. T. Miyagawa et. al. Physics-informed Neural Networks for Functional Differential Equations: Cylindrical Approximation and Its Convergence Guarantees, https://arxiv.org/abs/2410.18153\n\nLooking forward to any future discussions if you have any questions/concerns."}, {"review_id": "lQsirGjYbh", "reviewer": "Reviewer_fiA5", "summary": "This paper proposes a new architecture for learning neural functionals inspired by the Riesz representation theorem and applies this architecture to neural surrogate modeling of Hamiltonian PDE systems. The authors show that their method outperforms baseline models such as FNO and U-Net on a variety of problems, including 1d and 2d PDEs.", "questions": "1. Did the authors perform ablations to justify their choice of the SIREN architecture?\n2. To better understand the timescale of solving these PDEs, how fast does it take numerical solvers to generate the same amount of data as in the inference step (Table 6)?", "limitations": "The authors discuss the limitations of their work in Section 5.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\nThe paper is well-written and well-motivated. The background theory is explained clearly for all audience members. The experiments follow a logical pattern, going from toy examples to 1d problems to 2d problems. The evaluation criteria goes beyond simple MSE and includes plots of the Hamiltonian over time.\n\nThe authors write that to the best of their knowledge, “this is the first time an architecture has been developed to learn function to scalar mappings.” I am only aware of one other prior work that has done this. The discriminator module in [1] uses a “neural functional layer” to map a function to a real number. \n\nDespite some overlap between this work and [1], the authors clearly put effort into identifying the important architectural considerations for the neural functional, and the application to Hamiltonian systems is novel to my knowledge.\n\n**Weaknesses:**\nThe authors should discuss similarities and differences between their method and the “neural functional” introduced in [1].\n\nWhile the experiments are thorough and detailed, the number of baselines could be expanded, particularly for the 2d shallow water equations where the difference in error between HNF, FNO, and U-Net is less than for the other problems. There are a variety of architectures that have been introduced in recent years that have been competitive with (or outperformed) FNO and U-Net in some problems.\n\nSince HNF is agnostic to the resolution, I would suggest adding experiments using different resolutions and grids. For instance, zero-shot super-resolution can be compared with FNO and other neural operators. Training on multiple grids or resolutions can also be performed. Out-of-domain generalization experiments to different grids would be interesting as well.\n\n**Minor notes:**\n1. Typo “its functional derivative can obtained” in line 38.\n2. Typo “IN” in line 144.\n3. In Figures 2 and 3, the red and pink lines may appear too similar for some readers.\n\n**References:**\n[1] “Generative Adversarial Neural Operators”", "was_revised": true, "final_justification": "My initial concerns were an expanded comparison with prior works, a comparison with more baseline methods, and experiments in super-resolution and changing grids at test time. The authors addressed my points thoroughly by providing experiments on baselines, discretization invariance, and numerical solver timing. For this reason, I decided to increase my score from my initial response.", "rebuttal": "Thank you for taking the time to write a review and for your suggestions. We have implemented many of these below:\n\n### Related Works\nThank you for providing the reference. Since submitting the paper, we have also found more works that introduce a neural functional parameterized by a kernel integral [1, 2] as well as related works from other reviewers [3, 4] that are similar. Since we arrived upon the idea independently, much of the language in our paper unfortunately suggests that it is a novel idea. We will work to remedy this by changing the title and framing of the work, as well as amending claims of discovering this. \n\nDespite the existence of prior works, there are some innovations that we make from independently developing functionals for Hamiltonian systems. Firstly, we make an effort to verify that kernel integrals can approximate functional derivatives with analytically understood systems. Secondly, we use neural functionals to do a PDE prediction task, which is enabled by Hamiltonian mechanics. Lastly, we investigate implementing the kernel as a conditioned neural field with more advanced architectures (FiLM/SIREN), which is empirically necessary for good performance (see ablation studies below). Specific to the referenced work from Rahman et. al., our approach introduces more complexity into the kernel, since they only work with a linear, local kernel parameterized by an MLP, and the referenced work does not make use of functional derivatives. \n\n### Baselines\nWe have included two additional baselines for 2D SWE, a Transolver [5] and PINO [6] baseline. We decided to compare against Transolver, since it introduces a new architecture (attention-based) and has shown empirically good performance, and PINO, to compare against an established, physics-informed method. The modified results for 2D SWE are:\n\n**Table 1: Rollout Error of Models on 2D SWE**\n| Model             \t| Sines | Pulse |\n|-------------------------|----------|-----------|\n| Transolver (3.6M) | 0.084$\\pm$0.0081 | 0.122$\\pm$0.0074 |\n| FNO (7M) | 0.057$\\pm$0.002 | 0.117$\\pm$0.0009 |\n| PINO (7M) | 0.053$\\pm$0.005 | 0.114$\\pm$0.0003 |\n| Unet (3M) | 0.010$\\pm$0.0014 | 0.042$\\pm$0.0006 | \n| HNF (3M) | 0.026$\\pm$0.0003 | 0.021$\\pm$0.0015 |\n\nWe find that adding a physics-informed loss to FNO (i.e., PINO), helps to improve the performance slightly. The loss formulation we employed was: $L = L_{data} + w*L_{PDE}$, where $L_{PDE}$ is the residual of the 2D SWE equations, calculated with finite differences. For a fair comparison, we do not perform test-time finetuning of PINO. We report results with $w = 0.0001$, which was the best-performing model after a hyperparameter sweep. \n\nSurprisingly, Transolver does not perform as well as other methods. We investigated this through an extensive hyperparameter sweep and different model sizes, but was unable to improve its performance beyond what is reported. As a whole, the results are consistent with prior works that suggest Unets, while simple, are a well-performing architecture [7] on regular grid problems (0.01 relative L2 error is also quite small, and is qualitatively indistinguishable from the ground truth).  \n\n### Discretization Invariance\nAs suggested, we experiment with querying models at different discretizations than its training resolution. The results are below, with the training resolution in bold. Rollout errors are averaged over the validation set, where validation labels are either downsampled or re-solved at a higher resolution using the same initial conditions. Additionally, for 1D Advection, we run an experiment where the model is queried on an unseen grid (x=[0, 32], n_x=256), after being trained on a grid (x=[0, 16], n_x=128), which requires generating new validation samples.\n\n**Table 2: Rollout Error of Models at Different Resolutions, 1D Adv**\n| Resolution ($n_x$) | HNF | FNO (du_dt) | Unet (du_dt)\n|--------------|--------------|--------------|--------------|\n| $64$ | 0.0029 | 0.0639 | 1.1882 | \n| $\\mathbf{128}$ | 0.0033 | 0.0639 | 0.0753 |\n| $256$ | 0.0055 | 0.0497 | 1.3177 |\n| $256, x=[0, 32]$ | 0.0021 | 1.3438 | 0.0723 |\n\n**Table 3: Rollout Error of Models at Different Resolutions, 2D SWE (Sines)**\n| Resolution ($n_x,n_y$) | HNF | FNO | Unet |\n|--------------|--------------| --------------| --------------|\n| $(64, 64)$ | 0.028 | 0.066 | 0.098 | \n| $\\mathbf{(128, 128)}$ | 0.029  | 0.066 | 0.008 | \n| $(256, 256)$ | 0.026 | 0.066 | 0.110 | \n\n**Table 4: Rollout Error at Different Resolutions, 2D SWE (Pulse)**\n| Resolution ($n_x,n_y$) | HNF | FNO | Unet |\n|--------------|--------------|--------------|--------------|\n| $(64, 64)$ | 0.035 | 0.173 | 0.173 | 0.145 |\n| $\\mathbf{(128, 128)}$ | 0.024 | 0.174 | 0.117 |\n| $(256, 256)$ | 0.019 | 0.174 | 0.225 |\n\nIn all cases, HNFs are capable of zero-shot super-resolution, as is possible with Neural Operators. Additionally, the error is approximately constant across discretizations and on a new, unseen grid for 1D Advection. Predictably, FNO models have nearly constant error across discretizations. Interestingly, for an extrapolated grid in 1D Advection, the model struggles to make predictions. For Unet models, the error increases when the grid spacing changes, however for a constant spacing $dx$, it is able to extrapolate to unseen grids (i.e. [0, 16] -> [0, 32] with $dx$ unchanged). \n \n### Ablations\nWe perform ablations to study the effect of different model choices on HNF performance, across all PDEs (1D Adv, 1D KdV, 2D SWE-Sines). We consider the effect of different kernels (linear/nonlinear), architectures (MLP/SIREN), conditioning mechanisms (concat/FiLM) and receptive fields (local/global). Since this design space is large, we consider the effects of incrementally adding model features, where the base model is specified as (Linear, MLP, concat, local), and each successful model is cumulative:\n\n**Table 5: Ablation Studies for Different PDEs**\n| Model | 1D Adv | 1D KdV | 2D SWE-Sines | \n|-------------------------|----------|-----------|-----------|\n|      | Rollout Err. | Corr. Time | Rollout Err.| \n| Base | 1.248 | 51.75 | 0.113 |\n| +Nonlinear | 0.028 | 73.00 | 0.034 |\n| +FiLM | 0.016 | 65.25 | 0.030 |\n| +SIREN | 0.003 | 73.25 | 0.026 |\n| +Global | 0.057 | 151.75 | 0.024 |\n\nIn general, using a SIREN vs. MLP architecture improves performance. We can see that using a nonlinear kernel is also important, as all Hamiltonians tested are nonlinear. In some cases, design choices such as using a global kernel or FiLM conditioning are also important. \n\n### Numerical Solver Timing\nWe update Table 6 with the runtime of a single step of numerical solvers, compared to neural solvers. Numerical solvers are run on a AMD Ryzen Threadripper PRO 5975WX 32-Core CPU due to lack of GPU compatibility. In the case of 1D Advection, the analytical solution is used. Furthermore, randomly sampled initial conditions can cause variability in the stiffness of the solution, therefore, solver runtimes are averaged across 10 samples. We note that the 2D SWE solver (PyClaw) uses a Fortran compiler and is optimized for hyperbolic PDEs, while the KdV solver is purely Python-based. Despite the additional overhead of backprop in HNFs, it is still modestly faster than a dedicated numerical solver. For a fairer comparison, the numerical solver would need to be coarsened to match the accuracy of the neural solvers (which are pretty accurate), however, this still serves as a point of reference and gives an understanding of the general timescales involved.\n\n**Table 8: Time/Step, in ms**\n| Model | Adv | KdV | SWE-Sines |\n|-------------------------|----------|-----------|-----------|\n| FNO | 0.083 | 0.091 |  0.967 | \n| Unet | 0.138  | 0.146 |  1.345 | \n| HNF | 0.126 | 0.228 |  4.547 | \n| Numerical Solver| 0.000 | 110.9 |  33.26 | \n\n### References\n1. J. Bunker et. al. Autoencoders in Function Space, https://arxiv.org/abs/2408.01362\n2. D. Huang et. al. An operator learning perspective on parameter-to-observable maps, https://arxiv.org/abs/2402.06031\n3. Z. Hu et. al. Neural Integral Functionals, https://openreview.net/pdf?id=aQuqw6eVKP\n4. T. Miyagawa et. al. Physics-informed Neural Networks for Functional Differential Equations: Cylindrical Approximation and Its Convergence Guarantees, https://arxiv.org/abs/2410.18153\n5. H. Wu et. al. Transolver: A Fast Transformer Solver for PDEs on General Geometries, https://arxiv.org/abs/2402.02366\n6. Z. Li et. al. Physics-Informed Neural Operator for Learning Partial Differential Equations, https://arxiv.org/abs/2111.03794\n7. G. Kohl et. al. Benchmarking Autoregressive Conditional Diffusion Models for Turbulent Flow Simulation, https://arxiv.org/abs/2309.01745v3\n\nIf there are any further questions/concerns, we would be happy to address them. Thanks!"}, {"review_id": "q65kjPS3uo", "reviewer": "Reviewer_snot", "summary": "This work introduces a novel neural network called a Neural Functional (NF). Rather than learning functions from inputs to outputs, it learns to map input functions to a single number — a functional. The authors apply this concept to construct a particular model called the Hamiltonian Neural Functional (HNF), which is capable of learning physical systems' energy-like values (Hamiltonians) that obey particular rules (PDEs in Hamilton form). Of particular importance is predicting how these systems evolve in time while keeping their total energy constant. The model exploits automatic differentiation to compute the evolution of the system. It is experimented with in simple examples, one-dimensional PDEs such as Advection and KdV equations, and a two-dimensional shallow water model. HNF is compared to various models such as MLPs, FNOs, and U-Nets in this work.", "questions": "### **Questions for Authors**\nIn addition to the weaknesses noted above, we ask the authors to consider the following questions:\n\n1. Can you explain when and why your model can learn complicated (nonlinear) functionals? Do you have any theory to back this up?  \n2. Why should we believe that using autograd on a neural network gives the correct kind of derivative (the functional derivative) in your setting?  \n3. Your method is built around learning Hamiltonians and using them to evolve systems over time while conserving energy. However, many real-world physical systems are not purely Hamiltonian — they may include friction, damping, external forcing, or other types of energy gain/loss. Can your framework be adapted to handle these types of systems? If not, how much does this limit the range of problems your method can be applied to in practice? It would be helpful to understand whether your approach can be extended to include non-conservative dynamics, or if it's restricted to idealized, energy-conserving cases only.  \n4. How does the choice of grid or resolution affect your results? Have you tested what happens when the inputs are more coarsely sampled?  \n5. Why didn’t you compare your model to other physics-based methods like DeepONet, PINA, or Symplectic Neural Operators?  \n6. What is the actual runtime, memory cost and training data demand of your model during prediction compared to FNO or U-Net, especially on the 2D shallow water test?  \n7. How do we know how much training data is enough?  \n8. How does one represent the impacts of distributed inputs like bed topography or diffusion coefficient in the equations?\n\nI'm putting Borderline Accept for now due to the novel idea and giving the authors the benefit of the doubt, but if the questions got no good answers I reserve the right to mark it down.", "limitations": "N/A", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 4, "strengths_and_weaknesses": "### **Strengths**\n- The paper explores a fresh idea — learning energy-like quantities (called Hamiltonians) directly as functionals, which is different from how most models are designed.  \n- The authors try to base their method on solid math (using something called the Riesz theorem), which helps explain how it works, at least for simple cases.  \n- The experiments show that the proposed model (HNF) does a good job keeping the system’s energy stable over time, which is important for long simulations.  \n- The authors test their model on various examples — simple mathematical expressions, one-dimensional equations, and a more complex two-dimensional fluid system, which gives a broad look at how well the method works.  \n\n### **Weaknesses**\n- The theory the paper relies on (Riesz theorem) only works for simple cases, but the authors use it for complex, nonlinear situations without clear justification.  \n- There’s no formal proof that the model works well or converges, not even for simple cases.  \n- Some similar past work is not discussed or compared to, including other models that also learn energy functions or solve physics-based problems.  \n- In some tests, like the 2D shallow water case, the model keeps energy well, but doesn’t predict the actual system behavior better than others like U-Net. Improvements or explanations are needed. \n- The model may be slow during prediction because it needs backward passes every step, but this isn’t clearly shown or discussed in the main results.  \n- The method requires the knowledge of the J operator, which was said to be trivial in the paper but may not be that trivial if we learn from real data and the real system is unclear. It also requires the preparation of derivative as training data `∂u/∂t` for training. The computing cost of this was not discussed.  \n- The impact of grid resolution (training and test accuracy, and compute time) was not discussed. For high-resolution simulations, this approximation may not be trivial. And we need to ensure the cost is not higher than the direct numerical simulation.  \n- The authors did not discuss if their system works for problems with parameters, i.e., shallow water equation with bed topography, or advection-diffusion equation with varying diffusion coefficient.", "was_revised": true, "final_justification": "Questions have been addressed thoroughly. Idea is novel. Results look promising.", "rebuttal": "Thank you for the review. Apologies for the brevity of tone, since there were many suggestions and questions:\n\n### Additional Theory\nWith the feedback of Reviewer a4TR, we have found additional insight that connects our approach with an approximation of functional derivatives.Consider functionals $F[u]: H \\rightarrow \\mathbb{R}$, where functions $u(x) \\in H$ are in a Hilbert space $H$ with inner product $\\langle\\cdot, \\cdot\\rangle_H$. $u(x)$ can be written as $u(x) = \\sum_{k=1}^\\infty \\langle u,\\phi_k\\rangle\\phi_k(x)$, with a basis $\\{\\phi_1, \\phi_2, \\ldots\\}$. Truncating the sum to a finite number of terms $m$ allows the cylindrical approximation of functionals [1]:\n\n$$\n\tf(a_1, a_2, \\ldots, a_{m}) := F[\\sum_{k=1}^{m}a_k\\phi_k(x))]\n$$\n\nwhere $a_k$ are basis coefficients $\\langle u,\\phi_k\\rangle$. With the cylindrical approximation, functional derivatives can be written as:\n\n$$\n\t\\frac{\\delta F}{\\delta u} \\approx \\sum_{k=1}^{m} \\frac{\\partial f}{\\partial a_k}\\phi_k(x)\n$$\n\nwhich converges to the true derivative as $m\\rightarrow\\infty$. To discretize the derivative, we define a set of $n$ points $\\mathbf{x} = [x_1, x_2, \\ldots, x_{n}]$, which results in the discrete approximation:\n\n$$\n\t\\sum_{k=1}^{m} \\frac{\\partial f}{\\partial a_k}\\phi_k(x) \\approx \\sum_{k=1}^{m} \\frac{\\partial f}{\\partial a_k}\\phi_k(\\mathbf{x}) = \\begin{bmatrix}\n        \t\\sum \\frac{\\partial f}{\\partial a_k}\\phi_k(x_1) \\\\\\\\\n       \t\\sum\\frac{\\partial f}{\\partial a_k}\\phi_k(x_2)  \\\\\\\\\n       \t\\vdots \\\\\\\\\n       \t\\sum\\frac{\\partial f}{\\partial a_k}\\phi_k(x_{n})\n     \t\\end{bmatrix}\n$$\n\nWe are interested if the gradient of a kernel integral can represent this discrete object. We skip to the final result, which uses kernel functions of the form $\\kappa_\\theta(x, \\mathbf{u})$, where the kernel depends locally on the coordinate $x$ and globally on function values $\\mathbf{u} = [u_1, u_2, \\ldots, u_n]$. \n\nThe gradient of a nonlinear, global kernel integral can be written as:\n\n$$\n\t\\nabla_u \\int_\\Omega \\kappa_\\theta(x, \\mathbf{u})u(x) \\text{d}x \\approx \\mu \\nabla_u(\\sum_{i=1}^{n} \\kappa_\\theta(x_i, \\mathbf{u})u_i) = \\mu\\begin{bmatrix}\n    \t\\frac{\\partial}{\\partial u_1} (\\sum\\kappa_\\theta(x_i, \\mathbf{u})u_i)\\\\\\\\\n   \t\\frac{\\partial}{\\partial u_2} (\\sum\\kappa_\\theta(x_i, \\mathbf{u})u_i)\\\\\\\\\n   \t\\vdots \\\\\\\\\n   \t\\frac{\\partial}{\\partial u_{n}} (\\sum\\kappa_\\theta(x_i, \\mathbf{u})u_i)\n \t\\end{bmatrix}\n$$\n\nwith quadrature weight $\\mu$. The expression $\\frac{\\partial}{\\partial u_{j}}(\\sum\\kappa_\\theta(x_i, \\mathbf{u})u_i)$ can be expanded: \n\n$$\n\t\\frac{\\partial}{\\partial u_{j}}(\\sum_{i=1}^n\\kappa_\\theta(x_i, \\mathbf{u})u_i) =  u_1\\frac{\\partial}{\\partial u_{j}}\\kappa_\\theta(x_1, \\mathbf{u}) + \\ldots +  u_j\\frac{\\partial}{\\partial u_{j}}\\kappa_\\theta(x_j, \\mathbf{u}) + \\kappa_\\theta(x_j, \\mathbf{u}) + \\ldots +  u_n\\frac{\\partial}{\\partial u_{j}}\\kappa_\\theta(x_n, \\mathbf{u}) \n$$\n\nIn general, we can construct a kernel $\\kappa_\\theta$ such that $\\frac{\\partial}{\\partial u_{j}}\\kappa_\\theta(x_i, \\mathbf{u})$ is different for $i=\\{1, \\ldots, n\\}$. An example is provided based on FiLM:\n\n$$\n\t\\kappa_\\theta(x_i, \\mathbf{u}) = [W^u\\mathbf{u}]^i(W^xx_i +b^x) + [b^u]^i\n$$\n\nwhere $W^u \\in \\mathbb{R}^{n\\times n}, W^x \\in \\mathbb{R}^{1x1}, b^u \\in \\mathbb{R}^{n},b^x \\in \\mathbb{R}^1$ and $ [\\cdot]^i$ is an operation that takes the $i$-th row of a column vector. In this case, $\\frac{\\partial}{\\partial u_{j}}\\kappa_\\theta(x_i, \\mathbf{u})$ is equal to a function $d_{ij}(x_i) = W^q_{ij}W^xx_i$. With deeper kernels, functions $d_{ij}$ can become more expressive, and we can also absorb the basis coefficients into $d_{ij}$.\n\nThe $j$-th row of the gradient vector can be written as a sum of functions $\\sum_{i=1}^{n} d_{ij}(x_i)$, which reduces to the cylindrical approximation when letting $d_{ij}$ be arbitrary. Letting $n \\rightarrow \\infty$ increases the number of functions $d_{ij}$ and coordinates $x_i$ to increase the number of bases to approach the true derivative as well as convert from a discrete to continuous representation. The cylindrical approximation makes no assumption of linearity, therefore, using a global, nonlinear kernel integral can represent a functional derivative. The kernel function must have enough inputs; if it restricted to a linear kernel $\\kappa(x)$ or a local kernel $\\kappa(x, u(x))$, the same calculation only yields a single basis function. Intuitively, the kernel needs to depend on enough inputs so its gradient can be expressive enough. \n\n### Discussion of Prior Works\nOur related works include a few works from molecular dynamics (MD) that learn energy functionals and use their derivatives to update the simulation. Reviewer a4TR also suggested many works from the DFT field that use functional approximation, which we will include. \n\n### Accuracy vs. Conservation of Energy\nIn general, adding inductive biases to conserve energy can help accuracy, but is not strictly true. For example, learning the identity map conserves energy (since the state does not evolve), but is not accurate. While HNFs are not better than the Unet, the additional inductive biases help it to extrapolate to unseen conditions, which is the axis upon which we seek to improve. \n\n### Prediction Speeds\nWe discuss this in Section 5 and in Appendix B.2 of the paper and in Table 8 of our response to Reviewer fiA5, where we run timing experiments of different models and numerical solvers. The overhead of backprop increases the computational cost with respect to other models and scales with the number of query points, which, as you mention, is a limitation.\n\n### J Operator and du/dt Preparation\nOn real-world systems, prior works sometimes make an ansatz for contributions to the dynamics that conserve energies or satisfy certain properties (for example in climate [2]). If the Hamiltonian structure is not known, we can make a guess for its form based on similar, identifiable systems, and model other dynamics separately. In general, preparing du/dt for training has negligible overhead. It uses finite differences, which is a constant number of operations per sample, and du/dt can be precomputed and cached for training as well. \n\n### Impact of Grid Resolution\nWe perform a set of studies to compare the accuracy of models trained on a given resolution and queried at coarser or finer resolutions, in Tables 2-4 in our response to Reviewer fiA5. At test-time, HNFs can be used at different grids and resolutions without significant changes in accuracy. As for the computational cost, we compare the cost of our model to numerical solvers in Table 8 of our response to Reviewer fiA5. While the cost of our model does scale with the input resolution, the cost of numerical solvers also scales with increasing resolution. \n\n### Effect of Coefficients or Sources/Sinks\nIt is possible to use our work with varying coefficients and source terms with changes to the Hamiltonian structure. We provide examples below by introducing coefficients to 1D Advection as well as varying topography to 2D SWE.\n\nIn the presence of coefficients, the Hamiltonian for 1D Advection is:\n\n$$\n\t\\mathcal{H}[u] = \\int_\\Omega -\\frac{c}{2}(u(x))^2 dx\n$$\n\nwhich recovers the equation $-c\\frac{\\partial u}{\\partial x} = \\frac{\\text{d}u}{\\text{d}t}$ To accommodate varying coefficients, we can train HNFs by inputting coefficients to the model: $\\mathcal{H}_\\theta(u, x, c)$. In cases where the coefficients can be factored out (i.e., $\\frac{\\delta \\mathcal{H}[u, c]}{\\delta u} = c\\frac{\\delta \\mathcal{H}[u]}{\\delta u}$), once can simply scale the model output $\\nabla_u\\mathcal{H}(u, x)$ by $c$.\n\nSpatially-varying beds $b(x, y)$ can accelerate or decelerate flows, which acts as a source term on the momentum balance. The shallow water equations become:\n\n$$\n\t\\partial_th + \\nabla\\cdot (\\mathbf{v}h) = 0, \\qquad \\partial_t\\mathbf{v} + \\mathbf{v}\\cdot\\nabla\\mathbf{v}=-g\\nabla h -gh \\nabla b\n$$\n\n$h(x,y)$ now represents the height above $b(x,y)$. This source term changes the Hamiltonian to:\n\n$$\n\t\\mathcal{H}^b[\\mathbf{u}] = \\mathcal{H}^b[v_x, v_y,h] = \\int_\\Omega \\frac{1}{2}h(v_x^2 + v_y^2 ) + \\frac{1}{2}gh^2(1+b) \\text{d}A\n$$\n\nThe operator matrix $\\mathcal{J}$ is left unchanged. One can check that applying Hamilton's equations recovers the modified shallow-water equations. When comparing this derivation to one with constant bathymetry, one can show that:\n\n$$\n\t\\frac{\\delta \\mathcal{H}^b[u]}{\\delta u} = \\frac{\\delta\\mathcal{H}[u]}{\\delta u} + \\begin{bmatrix}\n   \t0\\\\\\\\\n   \t0\\\\\\\\\n   \tghb\n \t\\end{bmatrix}\n$$\n\nTherefore, source terms can be potentially captured by adding terms to model outputs. In PDEs where a decomposition like this is not possible, the source term will need to be added to the model input to produce the output $\\nabla_u\\mathcal{H}_\\theta(u, x, b)$.\n\n### Extensions to Non-Conservative Systems\nSimilar to [4], we may be able to model the dissipative dynamics jointly with the Hamiltonian with the Dissipation $[D(u, x), H(u, x)]$. For continuum systems, the theory governing this is not well-understood. We can potentially make an ansatz about how the derivative of the dissipation relates to the dynamics $du/dt$ but this is not well-studied. \n\n### Comparison to Other Models\nWe run more baselines for 2D SWE, in Table 1 of our response to Reviewer fiA5, including a physics-informed baseline.\n\n### Training Data\nTo determine how much training data is needed, we train HNFs on datasets downsampled by a factor f = ½, ¼ or ⅛ . In general, increasing the dataset size consistently improves model performance. \n\n**Table 1: HNF Performance vs. Training Set Size**\n| Dataset Size | Adv | KdV | SWE-Sines |SWE-Gaussian|\n|-------------------------|----------|-----------|-----------|-----------|\n| f=1/8 | 1.320 | 71.00 | 0.0293 | 0.0243\n| f=1/4 | 0.302 | 71.25 | 0.0274 | 0.0230 |\n| f=1/2 | 0.075 | 95.25 | 0.0254 | 0.0227 |\n| f=1 | 0.003 | 151.75 | 0.0249 | 0.0216 |\n\n### References\n1. https://arxiv.org/abs/2410.18153\n2. https://arxiv.org/abs/2404.10024\n3. https://arxiv.org/abs/2201.10085"}, {"review_id": "NqilbIaXVO", "reviewer": "Reviewer_a4TR", "summary": "Building upon the Riesz theorem (Theorem 3.1), integral kernel operator (5), and neural field (7), the authors develop an architecture that represents linear (and potentially non-linear) functionals. A training algorithm to adapt it to functional Hamiltonian problems is also proposed. The experimental results demonstrate that the proposed model outperforms baseline models on several PDE tasks, while it is a bit computationally demanding.", "questions": "- [Question] Could you elaborate on the following statement and provide an example?\n    \n    > ... functional approximation and learning with its gradients may find other uses, such as in molecular dynamics or design optimization.  (lines 11-12)\n    > \n    > Beyond predicting (future) PDEs (states), functionals arise in many physical systems, opening interesting opportunities to apply and improve neural functionals for other meaningful uses. (lineas 41-42)\n    \n- [Question] (lines 118-119) Why did the authors adopt FiLM?\n    \n- [Question] Could the proposed model be extended to handle functional differential equations that involve functional derivatives?\n    \n- [Question] How robust are the experimental results to the grid size? Do you have any recommendations or practical insights on this?", "limitations": "Discussed in Section 5.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n\n- (Technical Contribution) The combination of the integral kernel operator and neural field effectively solves functional Hamiltonian problems. Leveraging neural networks to simplify functional analysis is a promising direction, especially given recent advances in the numerical analysis of functionals.\n    \n- (Broader Impact) The underlying idea of the proposed model is simple and general and has the potential to be applied to a broader range of research areas, such as density functional theory (DFT) and functional differential equations.\n    \n- (Clarity) The paper is easy to follow and well-written. Background of the Hamiltonian mechanics, functional analysis, and neural field is provided, making the paper reader-friendly.\n    \n- (Reproducibility) The code is provided for reproducibility.\n    \n\n## Weaknesses\n\nI hope my comments below help improve the quality of the work.\n\n- The paper would benefit from a discussion of several related studies on neural functionals that are currently missing, which will highlight the authors’ contribution. To be clear, I am not suggesting the addition of extensive experimental comparisons with other baselines; rather, I would appreciate a discussion of the similarities and differences between the proposed method and prior work.\n    \n    - [Miyagawa & Yokota, 2024] Miyagawa and Yokota propose a neural architecture that can handle functional derivatives and solve functional PDEs. They also prove the convergence of their approximated solution converges to the true solution. They meticulously cite related studies and provide thorough discussions, which I think are helpful to highlight the contribution of the authors’ manuscript. I encourage the authors to compare their work with this study and discuss pros and cons.\n        \n    - [Zhou+, 2023a] [Zhou+, 2023b] Zhou et al. propose a neural functional that takes the weight parameters of DNNs as inputs and outputs vectors. This construction is particularly relevant to applications involving implicit neural representations (INRs). I encourage the authors to discuss the strengths and limitations of this neural functional in the context of their own work.\n        \n    - [Pederson+, 2022] [FIedler+, 2022] [Li+, 2023] Density functional theory (DFT) is an important application of functional theory. [Pederson+, 2022] and [FIedler+, 2022]  are review papers on DFT and machine learning. [Li+, 2023] is an example that apply deep learning to DFT. I encourage the authors to discuss a potential impact of the work on this field.\n        \n    - Please see also the references therein and add relevant papers to the manuscript.\n        \n- The mechanism by which functional derivatives can be efficiently computed via automatic differentiation appears to be a well-established approach; fundamentally, it stems from the discretization of spatial coordinates, where functional derivatives reduce to partial derivatives. This idea has long been employed in DFT, where functional derivatives are computed by discretizing spatial coordinates—reducing them to partial derivatives—and approximating them using neural networks (e.g., see [Meyer et al., 2020], [Costa et al., 2022], [Medvidović et al., 2025], and Appendix A in [Miyagawa & Yokota, 2024]). To my understanding, such approximations converge to the true functional derivative as the grid size decreases.\n    \n    - That said, the proposed method is specifically designed for solving PDEs, which is a perspective not typically addressed in DFT. I would recommend that the authors acknowledge this similarity and include a discussion. Relevant aspects for comparison include computational efficiency, convergence rate with respect to grid size, and solution accuracy.\n        \n- I would not recommend using the name “Neural Functional“ for the proposed architecture, as it has already used in several studies [Zhou+, 2023] and may be confusing (e.g., (liens 2-4) “In this work, we derive and propose a new architecture, the Neural Functional, which learns function to scalar mappings.” ([pdf](zotero://open-pdf/library/items/TMCVY4J3?page=1)) ).\n    \n- The following statements could come across as an overstatement:\n    \n    > (1) … no current neural architecture can theoretically model this infinite-dimensional to scalar mapping (=functional), let alone have meaningful functional derivatives. In this paper, we introduce the Neural Functional, an architecture that can provably represent functionals in the linear case… (lines 34-37)\n    > \n    > (2) … this is the first time an architecture has been developed to learn function to scalar (function-to-scalar) mappings. (lines 318-319)\n    \n    because (1) a neural architecture that can theoretically model functionals has proposed by [Miyagawa & Yokota, 2024], and (2) the neural operator can be used to map a function to a scalar. Strictly speaking, the model proposed by  [Miyagawa & Yokota, 2024] requires the number of basis functions $m$ and the neural network size $d$ to be large for accurate infinite-dimensional approximation; however, the authors’ neural functional also requires limit operation ($\\Delta x\\rightarrow0$ and $d \\rightarrow \\infty$ ) for provably accurate estimation of functionals (lines 102-158 & figure 1: approximating the integral and parameterizing the kernel).\n    \n    - (By the way, $m \\rightarrow \\infty$ in [Miyagawa & Yokota, 2024] corresponds to increasing spatial resolution and, at a high level, aligns with $\\Delta x\\rightarrow0$ in the authors’ framework, to my understanding.)\n        \n    - Again, the paper benefits from discussing the differences from [Miyagawa & Yokota, 2024].\n        \n- (Section 6: Related Works) I encourage the authors to clarify the advantage of the proposed model over other neural PDE surrogates [50-67].\n    \n- Error bars are missing in Tables 1, 4, 5, and 6.\n    \n- Error bars (standard deviations) in Tables 2 & 3 are not very informative, as they are based on only three trials. I would like to suggest increasing the number of runs, report all results, and/or use the box plot to provide a more comprehensive view of the variability..\n    \n- [Minor typo] (line 119) a -> an\n    \n- [Minor typo] (line 144) IN -> In\n    \n- [Minor comment] I recommend that the authors refrain from overusing semicolons, as it makes some sentences ambiguous.\n    \n\n## References\n\n- [Miyagawa & Yokota, 2024] Miyagawa, Taiki, and Takeru Yokota. \"Physics-informed Neural Networks for Functional Differential Equations: Cylindrical Approximation and Its Convergence Guarantees.\" Advances in Neural Information Processing Systems 37 (2024): 72274-72409.\n    \n- [Zhou+, 2023a] Zhou, Allan, et al. \"Permutation equivariant neural functionals.\" Advances in neural information processing systems 36 (2023): 24966-24992.\n    \n- [Zhou+, 2023b] Zhou, Allan, et al. \"Neural functional transformers.\" Advances in neural information processing systems 36 (2023): 77485-77502.\n    \n- [Pederson+, 2022] Pederson, Ryan, Bhupalee Kalita, and Kieron Burke. \"Machine learning and density functional theory.\" Nature Reviews Physics 4.6 (2022): 357-358.\n    \n- [Fiedler+, 2022] Fiedler, Lenz, et al. \"Deep dive into machine learning density functional theory for materials science and chemistry.\" Physical Review Materials 6.4 (2022): 040301.\n    \n- [Li+, 2023] Li, Tianbo, et al. \"D4FT: A deep learning approach to Kohn-Sham density functional theory.\" ICLR 2023 (Top25%).\n    \n- [Meyer+, 2020] Meyer, Ralf, Manuel Weichselbaum, and Andreas W. Hauser. \"Machine learning approaches toward orbital-free density functional theory: Simultaneous training on the kinetic energy density functional and its functional derivative.\" Journal of chemical theory and computation 16.9 (2020): 5685-5694.\n    \n- [Costa+, 2022] Costa, Emanuele, et al. \"Deep-learning density functionals for gradient descent optimization.\" Physical Review E 106.4 (2022): 045309.\n    \n- [Medvidović+, 2025] Medvidović, Matija, et al. \"Neural network distillation of orbital dependent density functional theory.\" Physical Review Research 7.2 (2025): 023113.\n    \n\n## Review Summary\n\nThe paper is well-written, and the main idea is clearly presented. The applicability of the proposed method to other domains is noteworthy. While the authors acknowledge the limitations of their theory for nonlinear functionals, addressing this would likely require a separate study, given the substantial challenges inherent in nonlinear functional analysis. Although the computational cost of the proposed model is higher than that of the baselines, it appears acceptable.\n\nThat said, there remain critical issues regarding the comparison with prior work, which may affect the central claims of the paper. For this reason, I am currently inclined to recommend rejection. However, I remain open to revisiting my evaluation and look forward to further discussion with the authors.", "was_revised": true, "final_justification": "The paper is well-written, and the main idea is clearly presented. The applicability of the proposed method to other domains is noteworthy. While the authors acknowledge the limitations of their theory for nonlinear functionals, addressing this would likely require a separate study, given the substantial challenges inherent in nonlinear functional analysis. Although the computational cost of the proposed model is higher than that of the baselines, it appears acceptable.\n\nIn the original submission, all reviewers raised concerns that several important prior studies were not cited.\nThrough the rebuttal, the authors clarified their contributions in relation to prior studies: Although this paper has an overlap in its scope with several prior studies, the authors develop their method on the Hamiltonian formulation, which can be an independent contribution of this paper.\nFurthermore, the authors added error bars and increased the number of runs in their additional experimental results. Although some results do not appear statistically significant, the variances of the proposed method are sometimes smaller than baselines, which can be beneficial for practitioners. Additionally, the results help clarify the conditions under which the proposed method works well or not.\nI believe that the paper's quality has significantly improved through the discussions with the reviewers, and the paper's contributions have been much clarified.\n\nOn the other hand, I am still concerned about the statistical significance of the experimental results because several error bars (signal-to-noise ratios) are large, and some of the results, if not all, are inconclusive, potentially affecting reproducibility of the work. Also, while the authors effectively addressed many concerns raised by the reviewers, the paper requires a major revision due to the aforementioned lack of citations, as outlined in the responses for all reviewers.\n\nOverall, I think the contributions of the paper slightly outweigh the concerns, and I changed my score from 3 to 4 (Borderline accept).", "rebuttal": "Thank you for your detailed review and for the suggested related works. We discuss these and some additional points below:\n\n### Related Works\nThank you for providing such an extensive list of prior works. Many of these are new to us, so we appreciate the perspective and view they provide. Since the submission, we have also found other works that approximate functionals with an integral kernel [1,2], and other reviewers have also pointed out similar works [3, 4]. We arrived upon the idea independently, which unfortunately accounts for much of the language/framing used in the paper. We are happy to change this, especially since we have a novel perspective on using functionals in Hamiltonian mechanics to evolve PDEs. The related works will be expanded to include suggested prior works on approximating functionals and Section 3 will present intuition on their capabilities from a theoretical standpoint, rather than a derivation. Any language claiming to be the first work on approximating functionals will be changed. Beyond these general revisions, we also discuss some differences to specific works:\n\n#### Comparison to Miyagawa & Yokota, 2024\nMiyagawa and Yokota approximate scalar functionals by a SIREN that takes time and basis coefficients as input, and outputs a functional. In this form, it is only applicable to FDEs. A direct extension to PDEs without introducing Hamiltonian mechanics isn’t obvious, since functionals are usually not involved in solving PDEs. Ours (and prior works) also approximate functionals with an integral kernel to provide another perspective from Neural Operator literature and Riesz representations. Despite these differences, we are able to show that, with an appropriate choice of kernel function, the derivative of an integral kernel can be expressed as the cylindrical approximation of a functional derivative, which ties these two works together. The proof sketch is given above as a response to Reviewer snot. \n\n#### Comparison to Zhou et. al. and Naming of “Neural Functional”\nThe similarities to this family of works are likely in just the name “Neural Functional” and the general idea of inputting functions to networks. The authors are mostly interested in principled approaches to learn with the weights of other neural networks, for tasks such as predicting test-time accuracy of models or classifying INRs based on their weights. In our work, we mainly consider functionals as a mathematical object that maps between functions and scalars, and the applications of functional derivatives in PDE prediction. \n\nHowever, we agree that the overlap in terminology is confusing and will work to remedy this. Specifically, we will change the title of the work to “Hamiltonian Neural PDE Solvers through Functional Approximation” (or something similar) to emphasize the novel perspective of Hamiltonian mechanics, which allows functional approximators to solve PDEs. Since the contribution of our work is no longer centered around introducing neural functionals (the mathematical object), we don’t see a need to introduce a new term and will refer to them more generically as functional approximators, integral kernels, learning the Hamiltonian, etc. \n\n#### Comparison to Machine Learning for DFT works\nThese seem to be similar to the introduced method for approximating functionals, as they also employ a kernel integral (where the kernel is an energy density). Surprisingly, these works introduce kernel integration many years before they were developed as a theoretical foundation for neural operators. Beyond architecture, many works also make use of functional derivatives to approximate potentials for downstream use. We will make sure to mention these similarities, however our application of this paradigm to predict PDEs is still previously unexplored, due to the necessary introduction of Hamiltonian mechanics. Hamilton’s equations allow a mapping from scalar functionals to dynamical quantities of interest (du/dt) which differs from the usual Newtonian perspective in PDE works. Furthermore, learnable kernel integration in DFT arises since it resembles known equations and empirical functionals, however in PDEs, it is derived from operator approximation and Riesz representation theorems, which offers a perspective from functional analysis. \n\n### Computing Functional Derivatives with AD\nThank you for providing the relevant references. DFT is a field we are not familiar with, however, we recognize there are many similarities and will discuss them in our work. A similar line of work in molecular dynamics (MD) learns an energy functional and uses its spatial derivative (force) to advance a simulation in time. We made note of the MD domain in our related works, but will expand it to create a section about functional approximation/derivatives in other domains. \n\n### Overstatements\nThank you for pointing these out, we agree that the quoted lines are overstatements. We will remedy this by removing language that suggests that we introduce functional approximation, but rather use it as a tool to obtain neural PDE solvers that have better abilities to conserve Hamiltonians and extrapolate to new scenarios. Your intuition about the approach taken by Miyagawa & Yokota converging to ours is correct and is something we investigate in the above proof sketch and will discuss in the work. \n\n### Clarifying Relationship to Other Neural Surrogates\nThank you for the suggestion. We will expand on the benefits of our method (conserving energies, extrapolation) as well as some of the drawbacks (higher training/inference cost, more effort to derive Hamiltonian structure) with respect to common neural surrogates in the field. \n\n### Statistical Significance\nThank you for the suggestion. We have doubled the number of seeds (up to 6 seeds) in our experiments. Here is a comparison of one such table (Table 2 in the paper):\n\n**Table 1: Comparison of Mean/Std calculated across 3 or 6 seeds**\n| Model | Adv (3 seeds) | Adv (6 seeds) | KdV (3 seeds) | KdV (6 seeds) |\n|-------------------------|----------|-----------|-----------|-----------|\n| FNO | 0.81$\\pm$0.17 |  0.83$\\pm$0.17 |  69.5$\\pm$8.2 | 68.0$\\pm$10.5 |\n| Unet | 0.52$\\pm$0.35 | 0.40$\\pm$0.29 |  125.7$\\pm$8.5 | 134.0$\\pm$10.6 |\n| FNO (du/dt) | 0.044$\\pm$0.002 | 0.048$\\pm$0.007 |  75.5$\\pm$1.9 | 77.6$\\pm$3.6 |\n| Unet (du/dt) | 0.068$\\pm$0.029 | 0.057$\\pm$0.024 |  127.4$\\pm$4.9 | 113.3$\\pm$15.0 |\n| HNF | 0.0039$\\pm$0.0002 | 0.0039$\\pm$0.0008 |  150.9$\\pm$3.3 | 151.1$\\pm$3.0 |\n\nFor brevity we do not present all the updated tables here. While some numbers have changed, the overall findings remain the same. We also add these error bars into Tables 1, 4, 5, 6 (across 6 seeds). Here is one such updated table (Table 1 in the paper):\n\n**Table 2: Results for Toy Examples, with Mean/Std calculated across 6 seeds**\n| | Base | | | OOD | | | Disc | | |\n|-------------------------|----------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|-----------|\n| Metric | MLP | FNO | NF | MLP | FNO | NF | MLP | FNO | NF |\n| $\\mathcal{F}_l[u]$ | 2e-05$\\pm$2e-05 | 2.8e-04$\\pm$1.6e-04 | 3e-16$\\pm$2e-16 | 0.046$\\pm$0.074| 0.27$\\pm$0.24 | 2.13e-07$\\pm$4.21e-07 | 49.29$\\pm$11.05 | 49.70$\\pm$11.16| 0.74$\\pm$1.16|\n| $\\delta\\mathcal{F}_l/\\delta u$ | 0.083$\\pm$0.012| 0.066$\\pm$0.034 | 0.0012$\\pm$0.00053 | 0.11$\\pm$0.026 | 0.12$\\pm$0.077 | 0.0012$\\pm$0.00053 | 2.64$\\pm$0.029 | 2.70$\\pm$0.073 | 0.077 $\\pm$0.047 |\n| $\\mathcal{F}_{nl}[u]$ | 0.029$\\pm$0.025 | 0.016$\\pm$0.016 | 0.002$\\pm$0.0007| 6709$\\pm$1567 | 6493$\\pm$1614| 2908$\\pm$816 | 381$\\pm$256 | 322$\\pm$266 | 55$\\pm$50 |\n| $\\delta\\mathcal{F}_{nl}/\\delta u$ | 1.325$\\pm$0.388 | 2.097$\\pm$1.902 | 0.045$\\pm$0.041 | 1730$\\pm$288 | 1723$\\pm$295 | 1079$\\pm$231 | 92$\\pm$35 | 88$\\pm$37 | 35$\\pm$17 |\n\nFor chaotic systems (KdV) or nonlinear functionals, spread can be large, which inflates standard deviations.\n### Example about Future Applications\nIn design optimization, for example, we are trying to find a set of design variables (i.e. points that define an airfoil spline) that minimize or maximize some scalar functional (i.e., drag of an airfoil). Most design variables are a discrete representation of an underlying object. For example, an airfoil may be defined by 100, 200, …, 1000 spline points, so a direct map from the coordinates to the functional is usually infeasible or at least not generalizable to different discretizations. Parameterizing the functional as a kernel integral is a natural remedy to this problem, and the gradient calculated by AD could be potentially used to adjust initial design variables towards a more desirable state. \n\n### Adoption of FiLM\nWe perform a series of ablation studies to justify the use of FiLM, in Table 5 in the response to Reviewer fiA5. We find that it can improve the performance of HNFs in some problems, but in general this is an empirical choice. Other options for incorporating conditional information into $k(x, u)$ can be investigated in the future. \n\n### Extensions to FDEs\nThis is an interesting direction, and it may be possible. Replacing F[u] with a learnable kernel integral seems possible, but one potential challenge will be using two quadratures to approximate dF[u]/dt: one for integrating the kernel and another for integrating the functional derivatives. An interesting aspect of our approach is it can be used independently of a basis or basis coefficients. \n\n### Robustness to Grid Size\nWe find that HNFs are approximately discretization invariant, shown in Tables 2-4 in our response in Reviewer fiA5. After training on a given resolution, they can be queried at coarser or finer grids without significant changes in performance. \n\n### References\n1. https://arxiv.org/abs/2408.01362\n2. https://arxiv.org/abs/2402.06031\n3. https://openreview.net/pdf?id=aQuqw6eVKP\n4. https://arxiv.org/abs/2205.03017\n\nThank you for the review. If there are any additional concerns/questions, we would be happy to discuss them."}]

Author rebuttal:

Meta-review: This paper proposes a novel neural architecture grounded in functional analysis, combining the Riesz representation theorem, integral kernel operators, and neural fields to represent linear (and potentially nonlinear) functionals. A tailored training algorithm for functional Hamiltonian problems is also introduced. The technical contribution is significant, as it demonstrates how neural networks can be leveraged to simplify the treatment of complex functionals, which has been a longstanding challenge in numerical analysis.
The reviewers consistently praised the clarity of the paper. The exposition is well-structured, with sufficient background provided on Hamiltonian mechanics, functional analysis, and neural fields, making the paper accessible to a broad audience. The work also shows strong reproducibility, with code made available. Experimentally, the model achieves superior performance over established baselines on several PDE tasks, highlighting its effectiveness, though at a somewhat higher computational cost.
Beyond technical merit, the proposed approach has promising broader impacts. Its generality suggests applicability to diverse domains such as density functional theory and functional differential equations. Overall, the combination of novelty, rigor, clarity, and potential for broad influence outweighs the computational overhead.
Decision: Accept (poster)