Enabling Automatic Differentiation with Mollified Graph Neural Operators
Authors: Ryan Y. Lin
Julius Berner
Valentin Duruisseaux
David Pitt
Daniel Leibovici
Jean Kossaifi
Kamyar Azizzadenesheli
Anima Anandkumar
Venue: Accepted by TMLR ()

Abstract: Physics-informed neural operators offer a powerful framework for learning solution operators of partial differential equations (PDEs) by combining data and physics losses. However, these physics losses require the efficient and accurate computation of derivatives. Computing these derivatives remains challenging, with spectral and finite difference methods introducing approximation errors due to finite resolution. Here, we propose the mollified graph neural operator ($m$GNO), the first method to leverage automatic differentiation and compute exact gradients on arbitrary geometries. This enhancement enables efficient training on arbitrary point clouds and irregular grids with varying geometries while allowing the seamless evaluation of physics losses at randomly sampled points for improved generalization. For a PDE example on regular grids, $m$GNO paired with Autograd reduced the L2 relative data error by 20× compared to finite differences, suggesting it better captures the physics underlying the data. It can also solve PDEs on unstructured point clouds seamlessly, using physics losses only, at resolutions vastly lower than those needed for finite differences to be accurate enough. On these unstructured point clouds, $m$GNO leads to errors that are consistently 2 orders of magnitude lower than machine learning baselines (Meta-PDE, which accelerates PINNs) for comparable runtimes, and also delivers speedups from 1 to 3 orders of magnitude compared to the numerical solver for similar accuracy. $m$GNOs can also be used to solve inverse design and shape optimization problems on complex geometries.

Published in Transactions on Machine Learning Research (10/2025)
Enabling Automatic Differentiation with Mollified Graph Neural Operators
Ryan Y. Lin', Julius Berner" Valentin Duruisseaux', David Pitt', Daniel Leibovici' Jean Kossaifi", Kamyar Azizzadeneshell", Anima Anandkumar'
Caltech NVIDIA
Reviewed on OpenReview:
https: // openreview. net/ forum? id- CGORINFAGr
Abstract
Physics-informed neural operators offer powerful framework for learning solution operators of partial differential equations (PDEs) by combining data and physics losses. However, these physics losses require the efficient and accurate computation of derivatives. Computing these derivatives remains challenging, with spectral and finite difference methods introducing approximation errors due to finite resolution. Here, we propose the mollified graph neural operator (mGNO), the first method to leverage automatic differentiation and compute enact gradients on arbitrary geometries. This enhancement enables efficient training on arbitrary point clouds and irregular grids with varying geometries while allowing the seamless evaluation of physics losses at randomly sampled points for improved generalization. For a PDE example on regular grids, mGNO paired with Autograd reduced the L2 relative data error by 20x compared finite differences, suggesting it better captures the physics underlying the data. It can also solve PDEs on unstructured point clouds seamlessly, using physics losses only, resolutions vastly lower than those needed for finite differences to be accurate enough. On these unstructured point clouds, mGNO leads to errors that are consistently 2 orders of magnitude lower than machine learning baselines (Meta-PDE, which accelerates PINNs) for comparable runtimes, and also delivers speedups from 1 to 3 orders of magnitude compared to the numerical solver for similar accuracy. mGNOs can also be usod to solve inverse design and shape optimization problems on complex geometries.
1 Introduction
PDEs are critical for modeling physical phenomena relevant to scientific applications. Unfortunately, numerical solvers become very expensive computationally when used to simulate large-scale systems. To avoid these linitations, neural operators, machine learning paradigm, have been proposed to learn solution operators of
PDEs (Azizadennshelli et al. 2024: Berner et al. 2025 Neural operators learn mappings between function spaces rather than finiteddimensional vector spaces, and, in particular can be used to approximate solution operators of PDE families (Li et al. 2020b Kovachki et al. 2023). One example is the Fourier Neural Operator (FNO) (Li et al. 2020a), which relies on Fourier integral transforms with kernels parameterized by neural networks. Another example, the Graph Neural Operator (GNO) (Li ct al. 2020b), implements kernel integration with graph structures and is applicable to complex geometries and irregular grids. The GNO has been combined with FNOs in the -eofmtry-Informed Neural Operator (GINO) (Li et al.) 2023) to handle arbitrary geometries when solving PDEs. Neural operators have been successfully used to solve
PDE problems with significant speedups (Kurth et al. 2022 Gopakumar et al. 2023 Wen et al.] 2023 Ghafourpour et al. (0255). They have shown great promise, primarily due to their ability to receive input. functions at arbitrary discretizations and query output functions at arbitrary points, and also due to their universal operator approximation property (Kovachki al. 2021
Work done partially Caltech.

---

Published in Transactions on Machine Learning Research (10/2025)
Mollified GNO Kernel (1.024 points)
GNO Kernel (1,024 points)
1.0-
bump
quartic
octic
half. cos
0.8
0.6
weights
Mollified Kernel (256 points)
points)
0.4
0.2
0.25
0.0
|| vl
Figure 1: (left) Mollisied GNO kernel's neighborhood and weights versus those of the vanilla GNO kernel
with differing point densities. (right) Examples of weight functions (equation 9 equation 12 for mGNO.
However, purely data-driven approaches may underperform in situations with limited or low-resolution data (Li et al. 2021b), and may be supplemented using knowledge of physics laws (Karniadakis et al.
2021
as additional loss lerms, as done in Physics-Informer Neural Networks (PINNs) (Rassi el al.] |01 7a 2019
In the context of neural operators, the Informe- Neural Operator (PINO) (Li et al. 2021b combines training data (when available) with a PDE loss at a higher resolution, and can be finc-tuned on a given PDE instance using only the equation loss to provide a nearly zero PDE error al all resolutions.
A major challenge when using physics losses is to efficienlly compute derivatives without sacrificing accuracy, since numerical errors on the derivatives compounded in the physics losses. Approximate derivatives be computed using finile differences, but can require very high resolution grid to be sulficiently accurate, thus becoming intractable for fast-varying dynamics. Numerical derivatives can also be computed using Fourier differentiation, but this requires smoothness, uniform grids, and performs best when applied to periodic problems. In contrast, automatic differentiation computes exact derivatives using repeated applications of the chain rule and scales better to large-scale problems and fast-varying dynamics. Unlike numerical differentiation methods, which introduce errors, automatic differentiation gives exact gradients, ensuring the accuracy required for physics losses, making it the preferred approach for chysics-informee machine learning.
Approach. We propose fully differentiable modification to GNOs to allow for the use of automatic differentiation when computing derivatives and physics losses, which was previously prohibited by the GNO's diffeetferntii.blllty. More precisely, we replace the non-differentiable indicator function in the GNO kernel integration by a differentiable weighting function. This is inspired by mollifiers in functional analysis Evans 2010), which are used to approximate, regularize, or smooth functions.
As a result, our differentiable mollified GNO (mGNO) is the first method capable of computing exact derivatives at arbitrary query points. The resulting mGNO can, in particular, be used within GINO to learn elliciently and accurately the solution operator of families of large-scale PDEs with varying geometries without data using physics losses. A data loss can help further improve the results when reference simulations are available, even if these are of very low resolution, when paired with a high-resolution physics loss.
Summary of Results. We test our approach on Burgers' and Navier Stokes equations with regular grids, and on nonlinear Poisson and hyperelasticity equations with varying domain geometrics. Figure shows example solutions, highlighting the complexity of the geometries considered. Physics losses prove critical, emphasizing the need for accurate and efficient derivative computation.
2

---

Published in Transactions on Machine Learning Research (10/2025)
Burgers
Nonlinear Poisson
Hyperelasticity
Figure 2: Examples of solutions for the problems considered.
Using Autograd instead of finite differences reduces the relative L2 data loss by 20x for Burgers' equation on regular grids, showing that the Autograd physics loss captures the underlying physics more accurately, In hybrid training with noisy reference data, the Autograd- based mGNO FNO remains robust and accurate, maintaining low PDE residuals and data loss.
It also excels at learning the lid cavity flow example governed by the Navier Stokes equations. For this problem, access to data is essential, but a small mumber of randomly sampled points (with Autograd) at each iteration achieves excellent results when trained with a hybrid loss. Including a physics loss also critical, it ensures PDE fidelity and guides the network toward physically consistent solutions oven low resolution.
Autograd mGINO performs seamlessly unstructured point clouds, while finite difference derivatives are not sufficiently accurate at the same training resolution and would need at least 9x more points. Autograd mGINO achieves relative error 2-3 orders of magnitude lower than the machine lcarning basclines (Meta- PDE 'Qin et al. 2022) using LEAP and MAML) considered for a comparable running time, and enjoys speedups of 20-25x and 3000-4000x compared to the numerical solver for similar accuracy on the Poisson and hyperclasticity equations. Furthermore, as a result of its differentiability, mGINO can be used scamlessly for solving inverse design and shape optimization problems on complex geometries, as demonstrated with an airfoil design problem.
2 Background
2.1
Neural Operators
Neural operators compose linear integral operators X. with pointwise nonlinear activation functions a to approximate nonlinear operators. A neural operator is defined as
o K., + b.) + K. +z
(1)
where Q and P are the pointwise neural networks that encode (lift) the lower dimension function onto higher-dimenssonal space and project it back to the original space, respectively. The model stacks L layers of oW + KG + bx) where Wi are pointwise linear operators (matrices), Kr are integral kernel operators, by are bias terms, and are fixed activation functions.
The parameters of a neural operator consist of all the parameters in P. 2, W, Kr, by. Kossaifi et al. (2025) maintains comprehensive open-source library for learning neural operators in PyTorch, which serves as the foundation for implementation.
3

---

Published in Transactions on Machine Learning Research (10/2025)
Fourier Neural Operator (FNO). A FNO Li al. 2020a is neural operator using Fourier integral operator layers
(K())))(()==2) = Rs (Fuz)
(2)
where Rs is the Fourier transform of periodic function g parameterized by ø. When the discretizations of both the input and output functions are given on a uniform mesh, the Fourier transform F can be implemented using the fast Fourier transform. We note that the FNO architecture fully differentiable. A depiction of the FNO architecture is provided Figure
14
Graph Neural Operator (GNO). A GNO (Li et al. 2020b) implements kernel integration with graph structures and applicable complex geometries and irregular grids. The GNO kernel integration shares similarities with the pressag-passingg implementation of graph neural networks (GNNs) (Battaglia al.) 2016). However, GNO defines the graph connection in ball in continuons physical space, while GNN assumes fixed set of neighbors in discrete graph. The GNN nearest-neighhor connectivity violates discretization convergence and degencrates into pointwisc operator at high resolutions, leading to poor approximation of the operalor. In contrast, GNO adapts the graph based on points within a physical space, allowing for universal approximation of operators.
Specifically, the GNO acts on an input function v as follows,
(cno(U)(z) = S 1 In.co) dy,
(3)
where C Idd is the domain o, k a learnable kernel function, and (x) is the indicator function over the ball (x) of radius centered at a E D. The radius r is hyperparameter, and the integral can be approximated with Riemann sum, for instance.
Geometry-Informed Neural Operator (GINO). GTNO (Li et al. 2023) proposes combine an FNO with GNOs handle arbitrary geometries. More precisely, the input passed through three main neural operators,
GGINO
GNO
GeNO
gencoder. GGNO
(4)
First, a GNO encodes the input given on an arbitrary geometry into latent space with regular geometry. The encoded input can be concatenated with signed distance function evaluated on the same grid if available. Then, an FNO is used mapping that latent space for efficient global integration. Finally, GNO decodes the output of the FNO by projecting that lalent space represenlation to the output geometry. As a result, GINO can represent mapping from one complex geometry to another. A depiction of the GINO architecture is displayed in Figure 15
2.2 Physics- Informed Machine Learning
Physics- Informed Neural Network (PINN). In the context of solving PDEs, a PINN (Raissi al. 2017a b. 2019) is neural network representation of the solution of PDE, whose parameters are learned by minimizing the distance to the reference PDE solution and deviations from known physics laws such as conservation laws, symmetries and structural properties, and from the governing differential equations. PINNs minimize a composite loss
= Ldata + XC
5
where Ldata measures the error between data and model predictions while Cphysics penalizes deviations away from physics laws. PINN overcomes the need to choose a discretization grid of most nunerical solvers, but only learns the solution for a single PDE instance and does not generalize to other instances without further optimization. Many modified versions of PINNs have also been proposed and used solve PDEs numerous contexts (Jagtap et al. 2020 Cai et al. 2022: Yu et al. 2022). When data not available, we can try learn the solution by minimizing Lphysiis only.
4

---

Published in Transactions on Machine Learning Research (10/2025)
While physics losses can prove very useful, the resulting optimization task can be challenging and prone to numerical issues. The training loss typically has poor conditioning as it involves differential operators that can be ill-conditioned. In particular, Krishnapriyan et al. (2021) showed that the loss landscape becomes increasingly complex and harder to optimize as the physics loss cocfficient À increases. The model could also converge to a trivial or non-desired solution that satisfies the physics laws on the set of points where the physics loss is computed (Leiteritz & Pflüger 2021). There could also be conflicts between the multiple loss terms when their gradients point in opposite directions. Even without such conflicts, the losses can vary significantly in magnitude, leading to unbalanced backpropagated gradients during training (Wang et al., 2021), and thus to the contribution and reduction of certain losses being relatively negligible during training. Manually tuning the loss coefficients can be computationally expensive, especially as the number of loss terms increases. Note that strategies have been developed to adaptively update the loss coefficients and miligate this issue (Chen et al. 2018 Heydari et al., 2019 Bischof & Kraus 2021
Physics-Informer Neural Operator (PINO). In the physics-informed neural operator approach et al. 2021b), an FNO is trained with low-rcsolution simulations (or even without data) and high-rosolution physics losses, allowing for near-perfect approximations of PDE solution operators. To further improve accuracy at test time, trained models can be fine-tuned on given PDE instance using only the equation loss and provide near-zero error at all resolutions. PINO has been successfully applied to many PDEs (Song) et al. physical validity, and data elficiency in operator learning compared to purely data-driven methods.
2022 Meng et al. 2023 Rosoisky et al. 2023 In practice, the PDE loss vastly improves generalization,
Physics-informed approaches without physics losses. Various approaches enforce physics laws in surrogate models without using physics losses. This can be achieved, for instance, by using projection layers (Jiang et al. 2020 Duruisscaux et al.
2024 Harder et al. 2024), by finding optimal linear combinations
of learned basis functions that solve a PDE-constraines optimization problem (Negiar et al.) 2022 Chalapathi iet al. 2024), by leveraging known characterizations and properties of the solution operator as for divergence- free flows (Richter-Powel. et al. 2022 Mohan et al., 2023 Xing et al., 2024)), or Hamiltonian systems with their symplectic structure (Burby et al. 2020; Jin et al. 2020 Chen & Tao 2021 Duruisseaux et al., 2023 These approaches are particularly advantageous when the system is well-understoon and physical raint can be explicitly encoded, ensuring strict compliance with the underlying physics. They also tend to require less computational overhead than methods that rely on minimizing physics losses. However, such approaches are often tailored to specific PDE struct ures and are limited to dynamical systems with well-characterived methods offer greater flexibility and broader applicability across a
solutions. In contrast, physics-loss-basee wide range of PDEs, including systems whose solution properties are not fully understood. By incorporating physics- based loss terms during training, as in PINNs and PINOs, the governing equations serve as regularizers, offering greater flexibility, broad applicability across diverse PDEs, and the advantage of requiring only knowledge of the underlying equations. Hybrid strategies are also possible, in which certain constraints are enforced directly in the architecture while others are imposed through physics losses during training. For instance, in incompressible Navier-Stokes Stokes problems, physics loss can penalize deviations from the momentum equation, while the divergenco-free condition can be enforced explicitly in the model (c.g. using projections as done in Jiang et al. (2020) and Duruisseaux et al. (2024)).
Physics-Infused Transformer Architectures. Recent parallel research has explored embedding certain structures of PDEs into Transformer architectures. Transolver (Wu et al. 2024 Luo et al. 2025 introduces a physics-attention mechanism that adaptively groups mesh points with similar physical states, enabling the model to capture complex correlations across arbitrary geometries efficiently. Unisolver (Zhou et al. 2025 conditions Transformer on complete physics information, incorporating both domain-wise components, such as equation symbols and coefficients, and point-wise components such as boundary conditions, achieving broad applicability across diverse PDEs. PDEFormer (Ye et al. 2025) represents PDEs as computational graphs and combines symbolic and numerical information through a graph Transformer and implicit neural representation, allowing mesh-free solutions and efficient handling of varied PDE types. Collectively, these architectures illustrate an alternative approach to integrating physical knowledge into Transformer structures, offering com- plementary strategies for improving the accuracy and flexibility of machine learning surrogates for PDE solvers.
5

---

Published in Transactions on Machine Learning Research (10/2025)
2.3 Computing Derivatives
To use physics losses, a major technical challenge is to efficiently compute derivatives without sacrificing accuracy, since numerical errors on the derivatives will be amplified in the physics losses and output solution.
Finite Differences. A simple approach is to use numerical derivatives computed using finite differences. This differentiation method is fast and memory-effcient (it roquires O(n) computations for an n-point grid). However, it faces the same challenges as finite difference mumerical solvers: it requires fine-resolution grid to be accurate and therefore becomes intractable for multi-scale and fast-varying dynamics. On point clouds, the stencil coefficients in finite difference formulas vary from point to point and must be computed each time, as outlined in Appendix [E] adding to the computational cost. The errors in the resulting derivatives will also vary across the domain depending on the density of nearby points.
Fourier Differentiation. Fourier differentiation is also fast and memory-eficient to approximate derivatives as it requires On log n) given an n-point grid. However, just like spectral solvers, it requires smoothness, uniform grids, and performs best when applied to periodic problems. Fourier differentiation can be performed on non- non-uniform grids, but the computational cost grows to O(n2), and if the target function non-periodic or non-smooth, the Fourier differentiation is not accurate. To deal with this issue, the Fourier continuation method (Ganeshram et al. 2025 can be applied to embed the problem domain into larger periodic space, at the cost of higher computational and memory complexity.
Pointwise Differentiation with Autograd. Derivatives can be computed pointwise using automatic disferentiation by applying the chain rule to the sequence of operations in the model. Autograd (Maclaurin et al.) 2015) automates this process by constructing a computational graph during the forward pass and leveraging reverso-mode differentiation to compute gradients during the backward pass. Autograd is typically the preferred method for computing derivatives in PINNs for variety of reasons (Baydin et al. 2017
Unlike finite differences, which can introduce large discretization errors when not computed on fine enough grids, Autograd provides eract derivatives, cnsuring the accuracy that is critical in physics-informed machine learning, where derivative errors are amplisied in physics losses. By exact derivatives, we mean that automatic differentiation computes gradients precisely from the differentiable computation graph, accurate up machine precision. While the graph represents a discretized approximation of the continuous PDE, in practice, the computed loss and gradients are extremely accurate, and automatic differentiation remains highly reliable even when evaluated on a relatively small set of points, as demonstrated in our numerical experiments.
Autograd computes gradients in a single pass (while finite differences require multiple function evaluations). It provides exact derivatives regardless of the mesh resolution, making particularly advantageous for large-scale problems, where finite differences become computationally intractable.
Autograd can compute higher-order derivatives with minimal additional cost, while finite differences require additional function evaluations and can suffer from further error accumulation.
Autograd can compute derivatives at any point in the domain scamlessly, while finite differences can struggle to handle complex geometries.
Figure 3 empirically shows that Autograd derivatives are smoother and more stable than finite differences with the proposed model.
However, Autograd also has limitations: all operations need to be differentiable, and storing all intermediate compulations in a computational graph can significantly increase memory usage for deep models. As a result, it can be slower and more memory-intensie for which finite differences are accurate enough at low resolutions. When the physics loss involves deep composition of operations, issues of vanishing or exploding gradients can also be exacerbated during backward propagation with automatic differentiation.
than finite differences for simpler low-dimensional problems

---

Published in Transactions on Machine Learning Research (10/2025)
Groundtruth
Autograd
Fifite-Differnnces
100
100
100
100
1.0
1.0
0.0
0.5
0.5
1.0
Figure 3: Qualitative assessment of the smoothness of derivatives. We fit mGINO with physics loss to the differentiable function u(r, y) = sin(4 on [0, 1]2, and compare the analytic ground truth derivative Oznu(s, y) = yy sin(Arwy) (left) with the numerical derivatives obtained using automatic differentiation (middle) and finite differences (right).
3 Methodology
Mollified GNO. Recall from Equation (3) that the GNO computes the integral of an indicator function [n)( that is not differentiable. We propose fully differentiable layer that replaces the indicator function (x) with a differentiable weight function w supported within (r)x). This is inspired by the mollifier in functional analysis (Evans 2010), which is used to smooth out the indicator function. The resulting mollified graph newral operalor (mGNO) acts on an input function via
Lult, u(x ())(, v)v(y)
Gm))NOL)(x))
(6)
for x D. Here, we padded the input of the FNO such that its output function is supported the extended domain
D= D = fr+ x D, y (,))),
which allows the recovery of exact derivatives at the boundary. Then, we can compute the derivative
O₂ (mcno(t))s)= Oz(w(x, v)s(,, v)]u(y) dy.
B..))
(7)
Automatic differentiation algorithms can then be used compute the derivatives appearing in physics losses. We can also use a cached version of neighbor search (i.c. store neighbors with nonzero weight) keep the method efficient the original GNO up to the negligible cost of evaluating the weight function w. We emphasize that these definitions work for arbitrarily complex domains D. As for the GNO, the radius is hyperparameter and the integral can be approximated with Riemann sum, for instance. An example of simplified pseudocode for the mGNO layer is provided in Appendix D
Mollified GINO. The mollified GNO can then be used within a fully differentiable mollified geometry- informed neural operator (mGINO) to learn efficiently the solution operator of PDEs with varying geometries using physics losses where the derivatives are computed using automatic differentiation,
GmGINO = deccdder o GeNO o gencder
(8)
This allows mGINO to be used for solving inverse design and shape optimization problems on complex geometries, by backpropagating gradients with respect to control parameters through it.

---

Published in Transactions on Machine Learning Research (10/2025)
Weight Functions. The choice of weight function plays central role in mollifying the hard cutoff neighborhoods of standard GNOs with smooth, compactly supported functions that enable stable and accurate automatic differentiation. These weights gradually attenuate contributions near the neighborhood boundary, improving gradient stability and reducing artifacts. Key properties include smoothness to avoid spurious gradient oscillations, compact support to control computational cost. by ignoring distant points, and locality tuned through the radius r.
In practice, smooth bell-shaped or cosine-like like profiles often perform well and are well-suited. The radius r is typically chosen to obtain desired typical number of neighbors per query point which balances efficiency and information aggregation. Performance is poor when the radius is too small, and when becomes too large, the computational and memory costs siguificantly increase while performance deteriorates. In practice one can start from a smaller value and progressively increase it until satisfactory performance achieved.
Letting d = || yll2/r, examples of weight functions w with support in (x) are
Wbump (I, 9) (z) (v) exp(d2/(d) -1)) (qquatti((I,))) = (v) (1 2d + d'))
(9)
(10)
(Feetic(1,,1) 2) = =BB.(x) (1 6d + 8d 3d"), Whalf_ Whalf_oos(1,,]] (z) (v) [0.5 + 0.5 cos((T())))))
(11)
(12)
These weight functions, displayed on the right in Figure are decreasing functions from 1 to 0 on [0, and preserve differentiability for all points.
An ablation study for the radius r and choice of weight function w conducted for the nonlinear Poisson equation in Section 1.2.6] and a discussion on the choice of aggregation scheme for neighbors contributions is provided in Section 4.2.3 in particular to highlight how mean aggregation can lead undesired artifacts when training based on data loss only.
Remark. The entire mGNO and mGINO architectures are designed such that all of their components are differentiable, ensuring that derivatives of any order can be computed correctly and accurately using automatic differentiation. For instance, we avoid non-differentiable activation functions such as ReLU (Rectified Linear Unit), which introduce discontinuities in their derivatives, and instead use smooth alternatives like GeLU (Gaussian Error Linear Unit). All the other operations in the network, including the Fourier layers, pointwise linear neural networks, and pointwise linear operators, also implemented in a differentiable manner. This design eliminates issues related to subgradients or undefined derivatives at specific points and preserves the validity of higher-order derivative computations using automatic differentiation. A discussion of using Autograd with non-differentiable components is included in Appendix G
4 Numerical Experiments
We use Meta-PDE Qin et al. 2022) and the popular finite element method (FEM) FEniCS (Alnæs et al. 2015 Logg et al. 2011 baselines. Meta- PDE learns initializations for PINNs over multiple instances that be fine-tuned on any single instance, and two versions have been proposed based on the meta-learning algorithms MAML Fiin el al. 2017) and LEAP (Flennerhag et al. 2019). While FEniCS has been successfully applied in various disciplines, it has not been highly optimized for our applications and could possibly be outperformed by other FEMs and non-FEMs [Liul 2009) such as Radial Basis Function (RBF) methods (interpolation using RBFs), Finite Point Methods (FPMs) and Moving Least Squares (MLS) methods (weighted least squares approximate solutions) and spectral methods (expanding the solution in basis of functions).
8

---

Published in Transactions on Machine Learning Research (10/2025)
4.1 Burgers' Equation
4.1.1 Problem Description
The Burgers' equation models the propagation of shock waves and the effects of viscosity in fluid dynamics. We consider the 1D time-dependent Burgers' equation with periodic boundary conditions, and initial condition uo E (D; R) with D = 1). The goal is to learn the mapping G from the initial condition u(æ, 0) 40
to the solution u(z,t) t) of the following differential equation for x E D:
Oeu(a,t) + Oz (z, t)/2) va t), u(x, uz)
tE (0.1)
We focus on the dataset of Li et al. 2021a consisting of 800 instances of the Burgers' equation with viscosity coefficient μ = 0.01 and 128 x 26 resolution. Each sample in the dataset corresponds to different
initial condition ug drawn from the Gaussian process (0, 625(-A + +51) Examples initial conditions and solutions to this time-dependent Burgers' equation are shown in Figures and
0.4
0.2
u(3r 0)
0.0
-0.2
Figure 4: Examples of initial conditions used for the Burgers' Equation.
Figure 5: Ground truth solutions for several initial conditions of the Burgers' Equation
The input initial condition u²(x) given on a regular spatial grid is first duplicated along the temporal dimension to obtain 2D regular grid, and then passed through a 2D FNO and a mollified GNO produce predicted function
== (GmGNo GFno)(u)
(13)
approximating the solution u(æ, t). We minimize weighted sum of the PDE residual and initial condition loss,
(Rungers(v) =118,v + 0.(v212) - v + a (l)(:.0) wolliR(D)
(14)
9

---

Published in Transactions on Machine Learning Research (10/2025)
4.1.2 Comparison of Differentiation Methods and Training Strategies
We trained models using Lurggrss where the derivatives are computed using various differentiation methods. The results, obtained by fine-tuning the models one instance at a time and then averaging over the dataset, are presented Table | Compared to numerical differentiation, mGNO FNO with Autograd was easier lo tune and produced a model with data loss 20x lower (although 6x slower per training epoch), suggesting it better captures the physics underlying the data. It also achieved physics loss x000 smaller, but this expected since the physics loss is evaluated using Autograd.
Table 1: Data and physics losses (computed using Autograd) for the proposed mGNO o FNO models and PINO baselines, trained with Burgers where derivatives are computed using different methods.
PDE Residual Loss
Relative L2 Data Loss
Autograd (ours)
1.62 10-6
1.33 10-2
Finite Differences
4.89 10-3
2.24 1001
Fourier Diflerentiation
4.77 10-3
2.19 1001
PINO Baselines:
Autograd
1,87 10-9
1.22 10
Finite Differences
3.71 10
2.26 10
Fourier Differentiation
3.63 10
2.19 10-1
Transolver Baseline (fiite-diifoornnee))
3.01 - 10-2
7.47 10-2
A comparison to PINO (i.c. using an FNO instead of mGNO o FNO) shows that PINO achieved a lower physics loss, but failed to reduce the Relative L2 data loss below 10% on the hyperparameter sweep considered. Note that Autograd mGNO FNO was only 1.2x slower than Autograd PINO (although this was with batch size of 1 and Autograd PINO could be accelerated more easily using batching). We also compared performance to Transolver (wu et al. 2024 with finite differences. We performed a grid scarch over lcarning rate, hidden dimension, number of heads and layers, initial condition loss coefficient, and reported the best result in Table [ Transolver does not reduce the PDE residual as effectively as the other approaches, but strikes a better balance by achieving lower data loss than the remaining baselines. Nonetheless, it still falls short of the performance of Autograd mGNO FNO, which attains substantially lower PDE residual while staying closer to the reference data.
We trained m&NOO FNO models using data loss, physics loss, and a hybrid loss (data and physics). The results displayed in Table show that using the data loss helps reduce the gap with reference data, but can come at the cost of higher physics loss, and vice versa.
Table 2: Data and PDE loss of three models trained with different losses C, using Autograd to compute the physics losses, for the Burgers' equation.
Ldsta
Laata + *Cphysics
Cphysies
Relative L2 Data Loss
1.42 10
3.81 10 -3
1.33 10-2 2
PDE Residual Loss
4.80 10
1.48 10-4
1.62 10-6 6
We also tried pre-training the mGNO FNO models using the data loss before fine-tuning them with the physics loss Burgers only, but were not able to obtain better results this way. The data loss rapidly deteriorates from its original value obtained in data-driven way, while the physics loss improves slowly but does not get better than when training with hybrid and physics-only losses without data loss pre-training.
10

---

Published in Transactions on Machine Learning Research (10/2025)
4.1.3 Subsampling and Data Efficiency
We now consider randomly subsampling the points at which the PDE residuals are evaluated to consider the effect of reducing the resolution used for the physics loss on performance. Table TableBshhon shows that despite the random locations and reduction of the number of points used evaluate the derivatives, we maintain good rosults until the number of points becomes very low.
Table 3: PDE Residual and relative data losses for mGNO FNO models trained with the physics loss Bargers using Autograd derivatives evaluated at different numbers of randomly subsampled points from the original 128 x 26 regular grid with 3328 points.
PDE Residual Loss
Relative T,2 Data Loss
Full Grid (3328 points)
1.62 10-6
1.33 10
Random 1000 points
6.70 10-6
1.32 10
Random 500 points
8.15 10-5
1.33 10-
Random 250 points
1.44 10
Random 100 points
4.06-10-5
1.91 1.91-10-2
Random 50 points
2.19 10
Random 25 points
1.29-107
2.93 10-2
Random 10 points
3.49 10-4
4.62 10
We also investigate data efficiency of hybrid loss training, and more precisely how performance changes when varying the number of points at which the data loss and physics loss evaluated when training the proposed miNNO FNO models using hybrid training loss Ldata + XLphysics- Table shows the results obtained with mGNO FNO models on the time-dependent Burgers' equation. We see that performance remains unaffected until the spatial resolution falls below 16 or the temporal resolution falls below 5. This corresponds to 8x spatial subsampling and 5x time subsampling from the original 128 x 26 grid, and highlights the data efficiency of the proposed approach: instead requiring high resolution 128 26 simulations (3328 points) for the data loss, we can obtain results of similar accuracy by supplementing low resolution simulations (e.g. at 22 x 9, Le. 198 points) with a physics loss evaluated at higher resolution.
Table 4: Relative L2 data loss and physics loss (computed using Autograd) for the proposed mGNO FNO models on Burgers' equation, when training with a hybrid loss. We consider the effect of reducing the resolution at which the data loss is evaluated during training, starting from the original 128 26 resolution. Each row corresponds to different temporal resolution (subsampled from 26). We show in orange the settings where the data loss starts being affected, and in red the settings where the data loss is significantly worse.
different spatial resolution (subsampled from 128), while each column corresponds
Space
26
13
5
Time
Physies
Data
Physics
Data
Physics
Data
Physics
Data
128
1.48 10-1
3.81 10-3
1.85
10-1
4.69
10-3
1.77
10-1
4.80
10-3
1.11 -
6.84 1844-00-
64
1.73. 10-1
3.83
10-3
1.95. 10-1
4.89
10-3
1.86 10-1
4.82
10-3
1.17 10-4
-0-3
43
1.74 10-4
4.04
10-3
1.85 10
4.84
10-3
1.88
10-4
5002-10-3
1.-4-10-4
6.86-10-3
32
1.74 10-4
4.22
10-3
1.85.10-
4.89
10-3
1.82
10-4
-.91.10-3
1.17.10-4
6.84.10-3
26
1.82 10-4
4.21
1.77
10
4.82
1003
1.92
10
4.87
1.13.10-4
6.88
22
1.70 10
4.03-10
1.68
10
3 4.60
1.72
10
4.77
1.21
7.03
16
1.84.00 10
8.60-10
1.99 10
8.08
10
1.97
10
7.27-10
1.41
8.33-10
11
1.86 10
1.32-10-2
1.90 10
1.17-10
2.04
10
1-11-10-2
1.40
1.13-10-2
7
1.57
10
5.67
1.68 10
5.15-10-2
2.01
10
5.06
10-2
1.21-10
5.12-10-2
11

---

Published in Transactions on Machine Learning Research (10/2025)
4.1.4 Robustness to Noisy Data
We now investigate the effect of noisiness on reference data when training the proposed mGNO o FNO models using a hybrid training loss Latta + Xh both with finite differences and automatic differentiation.
More preciscly, instead of the distance to the reference solution u, the data loss is based on the distance to noisy reference solution
(. where 17 is the noise level, the integral represents the average absolute value of u(2, t), and e(x, t) is an independent standard normal random variable for each (r, (), that E i..d 1N.
ū(x, u(z,)) -nn--neeee e(x,t) S.-o lu(a' (f), dt',
(15)
Table 5 illustrates the effect of noise in the reference data on the mGNO FNO models trained with a hybrid physics-data loss. The results show that the Autograd approach, our proposed method, maintains very low PDE residuals even noise increases, demonstrating ils robustness and stability. The relative L2 error with respect to the true solution grows with noise level, reflecting the unavoidable discrepancy between noisy training data and the true solution, but increases significantly slower than the training L2 error with respect to the noisy reference. This indicates that the physics loss truly helps the model to learn effectively the dynamics without overfitting to noise.
Compared to finite differences, Autograd achieves slightly lower residuals and errors across noise levels, highlighting its reliability for accurately computing derivatives within the physics loss. These results confirm that the proposed Autograd-based hybrid training is effective and resilient, preserving model accuracy even when the reference data is noisy, while finite differences provide a reasonable but less precise baseline.
Table 5: Effect of injecting noise of different levels nj (as defined in equation 15 in the reference data when training a mGNO o FNO on Burgers' equation with a hybrid loss. We report the PDE residual and the relative L2 data loss both to the true reference solution u and to the noisy reference solution û. We consider both the cases where the derivatives are computed using Autograd and finite differences in the physics losses.
Noise Level 1
PDE Residual
Relative L2 (true)
Training Relative L2 (noisy)
Autograd
1) 0% (Reference)
3.43 10-1
1.68 10-3
1.68 10-3
17 0.1%
3.39 10-4
1.70 10-3
1.90 10-3
77 0.5%
3.08 10
1.99 10
4.62 10 10-3 3
17 = 1%
2.79 10
2.64 10-3 3
8.74 10
? 2%
3.01 10
4.09 10-3
1.72 10-2
1] 5%
4.93 10-1
8.41 10-3
4.29 10-2
10%
5.73 10
1.31 10-2
8.58 10
77 20%
9.79 10-
2.22 10-2
10-1
50%
1.67 10
4.39 2
3.97 10
Finite Differences
77 0% (Reference)
4.78 10
1.35 10-3
1.35 10-3
0.1%
4.88 10-1
1.40 10-3
1.63 10-3
17 0.5%
-000
2.24 10-3
4.73 10-3
17 = 1%
10
3.54 10-3
9.10 -0
17 2%
6.28 1
5.81 5181
1.78 10-2
1] 5%
8.23 10-4
1.12-10-1
4.40 10
10%
1.16 10-3
1.85 10-2
8.74 10-2
17 20%
1.71 10
3.07 2
1.72 10
50%
2.98 10-3
6.08 10-2
4.01 10
12

---

Published in Transactions on Machine Learning Research (10/2025)
4.2
Nonlinear Poisson Equation
4.2.1
Problem Description
The Poisson equation is a fundamental PDE that appears in numerous applications in science due to its ability model phenomena with spatially varying and nonlinear behaviors. We consider the nonlinear Poisson equation with varying source terms, boundary conditions, and geometric domain,
V [11 [11 0.1x()))))
fxx
E Ω
(16)
u(x) = b(x)
E an
(17)
where E R and Ω C R2, The domain Ω is centered at the origin and defined in polar coordinates with varying radius about the origin
(0)) ro 1 + e1 cos(40) + cos()0)),
where the parameters 01 and C2 are drawn from uniform distribution on (-0.2, 0.2). The source term f is of radial basis functions f(x) E 3r exp ||x p4l12, where B3 R and μ4 E R are both drawn from standard normal distributions. The boundary condition b is periodic function, defined in polar coordinates as bo + [b, cos(8) + by sff + by cos(20) + by sin() where the paramcters b; are drawn from a uniform distribution on 1), This is the setting used by Qin et al. 20022..
The mesh coordinates, signed distance functions, source terms, and boundary conditions, are passed through GINO model of the form
MGNN
GENO
GNO
Geneader
(18)
to produce approximation to the solution We minimize the PDE residual and boundary condition loss,
Proisson (v) = V (.1.) - -f + a [|l bll bl
(19)
Examples of predictions made by the trained mGINO model for variety of geometries are displayed in
Figure 6 below the corresponding reference solutions.
Ground Truth
Prediction
Figure Comparison of reference solutions the nonlinear Poisson equation on various geometries (top) with the corresponding mGINO predictions (bottom).
4.2.2 Comparison to Baselines
The results in Figure 7 show the trade-off between inference time and accuracy. mGINO achieves a relative squared crror 2-3 orders of magnitude lower than Meta-PDE (i.c. initialized PINNs) for a comparable running time, and a speedup of 20-25x compared the solver for similar relative accuracy. In addition, mGINO is more consistent across different instances, with smaller variations in errors compared to Meta-PDE.
13

---

Published in Transactions on Machine Learning Research (10/2025)
Nonlincar Poisson Equation
Hyperelasticit Equation
10a
NAMI
10P
LEAP
LEAP
10-1
mCINO [ours]
mCIINO [ouns
FExiCS
10-1
10-2
Delative Squared Eror
Relative Squard Eror
10-3
10-2
10-1
10-1
10-1
10-7
Time (3)
Time
(5)
Figure 7: Computational time for inference versus accuracy for the nonlinear Poisson and hyperelasticity equations. We compare the proposed approach with Meta-PDE (MAML and LEAP) and with the baseline FEniCS solver. The FEniCS solver was used with 6 different resolutions for the Poisson equation, and 4 different resolutions (itcratively trics to refine the solution) for the hyperclasticity equation. The mcan values across 200 instances are plotted, with the shaded regions representing the min and max error values.
4.2.3 Using Different Training Losses
We trained models using a data loss, physics loss, and a hybrid loss (data and physies). Table T shows that using a physics loss has a comparable data error to the data-driven approach, while achieving PDE residual 4-5 orders of magnitude lower. Fine-tuning the trained data-driven model using CPoisson did not work well, as its physics loss is 7 orders of magnitude higher than its data loss, and attempts at lowering the physics loss only worked by completely sacrificing the data loss.
Table 6: Data and PDE losses of mGINO trained with different losses & for the Poisson equation.
Ldata
Laata + phyyyscc
Ephysics
MSE Data Loss PDE Rosidual Loss
1.36 - 10-5
2.29 10
1.39 10
3.55 102
2.07 10-2
7.21 - 10-3
When training only with data loss, the predicted solutions have discontinuitics higher resolutions coinciding with circles of radius centered at the latent query points of the GINO (sce Figure 8 [a)) These discontinuities
lead to high derivatives and inaccurate physics losses, regardless of the differentiation method used (see Figures (c)). This happens despite the low data MSE, indicating that only training the mGLNO model with data loss is not sufficient to capture the solution correctly at higher resolutions.
Recall that the GNO's kernel integration can be viewed as an aggregation of messages we construct a graph on the spatial domain of the PDE, as described in Li et al. (2020b). The mean aggregation is given by
24 = a We((z) + [x)x
1
Σ
no(e(<z ())
vEN(()
where (z) R" are the node features, e(x,y) y) E R= are the edge features, W E Rax= is learnable, N(x) is the neighborhood of and ko(e(z, v)) is neural network mapping edge features to matrix in Ruxn When using a differentiable weight function, points in N(x)) at the edge of the neighborhood have ncar-zero weights but still contribute to the denominator IN(x)) Thus, as the query point x moves slightly, additional neighbors get included or excluded with near-zero weights, thereby introducing the discontinuitiss we see in Figure a). Using a sum aggregation for the output GNO's kernel integration mitigates the rigid patterns. When training with a physics loss, the patterns disappear, as shown in Figure g while the MSE for the predictions remains of the same order of magnitude.
14

---

Published in Transactions on Machine Learning Research (10/2025)
21
Predicted Solution
0.59
51
Gradient harm (uaing eutegrad)
01
Gredient Norm Joking diffemencee)
023
0.57
0
038
0
(a)
(b)
(c)
06
-21 50a
-01
-05
Figure (a) Prediction of mGINO trained using data loss only for the nonlinear Poisson equation. This is the predicted solution evaluated at high resolution small patch centered at latent query point of the mGNO. The prediction exhibits discontinuities that coincide with the circles of radius r (blue lines) centered at the neighboring latent query points. (b) Norm of the gradient of the predicted solution shown in (a), computed using automatic differentiation (in (b)) and finite differences (in (c)).
Predicted Solution
Ground Truth Solution
0.1
0.1
0.07
0.06
0.06
0.05
0.05
0.04
0.04
-011
0.1
Figure Comparison between the prediction of a mGINO model (left) trained using a hybrid loss, L Lasta phyysics and the corresponding ground truth solution (right) for the nonlinear Poisson equation. The prediction and ground truth solutions evaluated at high resolution small patch centered at latent query point of the mGNO. Using the physies loss regularized the higher resolution prediction, and removed the visible discontinuity patterns shown in Figure a).
4.2.4 Nonlinear Poisson PDE Residual Loss with Finite Differences
Figure 10 shows gradients norms domain patch, computed using Autograd and finite differences at 16 x 16 and higher resolutions (both displayed at 16 x 16). The finite difference gradient norms al lower resolution in (c) differ completely from those obtained using Autograd (b)(e), and high-resolution finite differences (f), while the latter are very similar. In addition, inaccuracies can be further amplified with finite differences on unstructured point clouds in regions with a low density of points. The Autograd mGINOs were trained successfully with a density randomly located points on the whole domain slightly lower than the density here with 162 points on the domain patch displayed in Figure 10 Given that numerical errors made derivatives amplified in physics losses, Lpoisson does not provide enough information to move towards physically plausible solutions when computed using finite differences unstructured point clouds at the same resolution at which Autograd mGINOs models were trained successfully. This showcases how Autograd m.GINOO can be used to oblain surrogale models using physics losses al resolutions al which finite differences are not accurate enough. We estimate that at least 9x more points would be necessary to compute reasonable finite difference derivatives. Higher resolution finite difference computations typically require a higher computational time and memory requirement (sce Appendix F for a comparison of training times different resolutions).
15

---

Published in Transactions on Machine Learning Research (10/2025)
Predicted Solution
Gradient Norm (using autograd on 16x16)
Gradient Norm (using FD on 16x16)
(a)
(b)
(c)
0.5
1.0
0.5
0.6
0.8
08
0.7
0.6
0.6
08
0.4
04
-0.
42
02
-1A
6.0
0.0
-0.4
Y
-0.4
-1.4
Predicted Solution
Gradient Norm (using autograd 96x96)
Gradient Norm (using FD 96x96)
(d)
(n)
-.05
-1.0
-35
Y
0.4
644
0.4
Figure 10: Prediction of mGINO model and the norm of its gradient, computed 16 x 16 regular grid (top row) and 96 x 96 regular grid (bottom row) using Autograd and finite differences (FD). All plots are displayed on a 16 x 16 grid for ease of qualitative comparison.
4.2.5 Data Efficiency
We now invostigate the data efficiency of hybrid loss training, and preciscly how performance changes when varying the number of points at which the data loss and physics loss are evaluated when using hybrid training loss Ldata + Xppysics for the nonlinear Poisson equation. Table 7 displays the results obtained with the proposed mGINO model, showing that both the data and physies losses worsen when the number of points used for the physics loss during training decreases. However, performance is not very heavily impacted when the number of points used for the data loss during training reduced.
Table 7: Data loss (MSE) and physics loss (computed using Autograd) for the proposed mGINO models on the Poisson equation, when training with hybrid loss. We vary the resolutions at which the data loss and physics losses are evaluated during training. Each row corresponds a different number of points used for the physics loss, while each column corresponds different number of points used for the data loss.
Data
400 points
200 points
50 points
Physics
Physics
Data
Physics
Data
Physics
Data
800 points
2.63 -0
8.28 10
3.31 10
3.89 10
3.82 10-2
6.73 10 10-5
400 points
2.13 10-1
1.10 10
1.22 10-1
9.20 10-5
1.39.
1.29 10-1
200 points
3.75 10
2.40 10
2.55 10-1
2.92 10-4
3.33 10-1
5.18 10-4
16

---

Published in Transactions on Machine Learning Research (10/2025)
4.2.6 Ablation Study on GNO Radius and Weight Function
We investigate how performance changes as the mGNO radius changes for mGINO models, and also compare the performance of the different weight functions (9) 12 on the nonlinear Poisson equation. Note that the radius cutoff used the weight functions of mGNOs the same as the mGNO radius to avoid the additional cost of carrying second neighbor search. The mGNO radius is an important hyperparameter to tune. It. needs to be large enough that the ball B, (z) contains sufficiently many other latent query points, and a value too large can lead to prohibitive computational and memory costs. We denote the number of latent query points in (x) by #IB,-xx)-- For the nonlinear Poisson equation, we are using a regular 2D latent space grid, so (x) will only contain a single latent query point (I itself) until r is large enough for Br(x)) to contain one extra latent query point in each direction. (x) will progressively contain more latent query points by thresholds as the radius increases further, together with computational time and memory cost.
Table shows the results obtained for mGINO models with different weight functions for values of such that (x),(x)]= 32,5², 72,99 15a. As expected, performance is poor when the radius is too small. On the other hand, when r becomes too large, the computational and memory costs significantly increase while performance deteriorates. For this numerical experiment, the best performance was achieved with radius for which (I), = =??, and Whalf cos and Equartic lod to the best performance.
Table 8: MSE of mGINO models trained on only PDE loss, with different GNO radii and weight functions.
bump
half cos
quartic
octic
(+)( 9
(r 0.0875)
3.64 10-1
3.41 10-1
2.20 10-2
3.48 10-1
(B) = 25
(r 0.13125)
4.51 10-5
3.59 10-5
3.95 10-5
4.24 10-3
# (z)(x) = 49
(r = 0.175)
1.40 10
1.39
2.93 2
9.49 10
#B, = -1
(r = 0.21875)
5.35 10-4
1.01 10-4
7.79
1.76 10-4
BB, = 225
(r = 0.355
3.92 10-3
9.36 10-4
4.28 10-3
6.59 10-4
4.3 Hyperelasticity Equation
4.3.1 Problem Description
The hyperclasticity cquation models the shape deformation under external forces of hyperclastic matcrials (e.g. rubber) for which the stress-strain relation is highly nonlinear. We consider the deformation of homogeneous and isotropic hyperelastic material when compressed uniaxially (assuming body and traction forces), and learn the final deformation displacement mapping the initial reference position to the deformed location. More precisely, we consider the deformation of a two-dimensional porous hyperelastic material under compression, as in Overvelde & Bertoldi (2014) and Qin et al. (202)). We keep the pores circular and fix the distance between the pore centers, so that the size of the pore is the only varied parameter. The size of the pores determines the porosity of the structure and affects the macroscopic deformation behavior of the structure. Examples of final deformation displacement fields are shown in Figure 11
Figure 11: Final deformation displacement for several instances of the hyperelasticity dataset.
17

---

Published in Transactions on Machine Learning Research (10/2025)
The mesh coordinates and signed distance functions passed through GINO model of the form
mGNN
GENO
Geneoder 9GNO
(20)
to produce approximation to the solution to the hyperclasticity equations.
The solution can be obtained the minimizer of the total Helmholtz free energy of the system Sa 1i dx. Here v denotes the Helmholtz free energy relating the Yola-KKrchhoof stress P with the deformation gradient de F via P= aF Instead of minimizing the PDE loss from the strong form of the hyperelasticity equation, we minimize the Helmholtz free energy of the system by randomly sampling collocation points from the PDE domain and the Dirichlet boundary and using these points form a Monte Carlo estimale of the total Helmholtz free energy. In addition, we add a weighted boundary loss term.
4.3.2 Numerical Results
We trained models using a data loss, a physics loss, and a hybrid loss. Table g shows that using a physics loss is critical for this problem, as the data-driven approach achieves a low data loss, but the predictions do not satisfy the physics at all. In contrast, both the physics-only and hybrid approaches achieve low data and PDE errors. As for the Poisson equation, fine-tuning the trained data-driven model using the physics loss did not work, due to the high physics loss of the data-driven model. As in Section 4.2.4 Autograd enables computation of accurate derivatives at resolutions where finite differences are not accurate.
Table 9: Data and PDE losses of mGINO trained with different losses & for the hyperclasticity equation.
Laata
Laata + AC
Cphysies
MSE Data Loss
3.35 10-7
9.69 10-7
9.42 10
PDE Residual Loss
2.97 1026
1.57 10-2 2
1.21 10
A comparison to baselines is displayed in Figure mGINO achicves a relative squared error orders of magnitude lower than Meta-PDE (i.e. initialized PINNs) with a slightly faster rumning time, and achieves a speedup of 3000-4000x compared to the FEniCS solver for similar accuracy. In addition, mGINO achieves consistent results across different samples (all within single order of magnitude) while Meta-PDE results span more than 3 orders of magnitude. This high variation in Meta-PDE results is likely caused by the difficulty disparity across samples, samples with larger pores are harder to resolve.
4.4
Navier-Stokes Equations
4.4.1 Problem Setting
Finally, we consider the lid cavity flow governed by the Navier Stokes equations
0u((2, t) u(x f) Vu(x, t) = Vp(x,t)+ Re Au(r,t),
E (0, .t. (0,
(21)
V . u(x t) = 0,
x (0, 1)2,) E [0,7]
u(z, 0) uo (2))
E (0, 1)²
We study the standard cavity flow on domain D (0, 1)2 with 10 seconds, where is velocity, p is pressure, P is density, and Re 500. We assume zero initial conditions, and the no-slip boundary condition where u(z, f) (0, 0) at the left, bottom, and right walls and u(I, t) = (1, on top. We also start from zero initial conditions, w(0 t) = 0 for x E 1)2. We are interested in learning the solution for the time interval [5, 10). Figure 12 shows the solution to this equation at a few timesteps. The main challenge lies in handling the boundary conditions within the Srbocitypreessro formulation.
18

---

Published in Transactions on Machine Learning Research (10/2025)
u velocity at = s
velocity at t= 7.5
velocity at t= 10.0
velocity 5.0
velocity tt 7.5
velocity at f 10.0
Figure 12: Example of solution for the lid cavity flow governed by the Savier-Stokee equations.
4.4.2 Results
We first trained mGNO FNO models using a data-only loss, a physics-only loss, and a hybrid loss. To facilitate the tuning of the loss coefficients for this numerical experiment, we employ ReLoBRaTLo (Relative Loss Balancing with Random Lookback) (Bischof & Kraus 2021) to adaptively update and balance the various loss coefficients in the total loss. ReLoBRaLo uses the history of loss decay across random lookback intervals to achieve uniform progress across the different terms and faster convergence overall. We also investigate how performance changes when varying the number of points at which the data and physics losses are evaluated during training when using a hybrid training loss Ldata + XCphysics for the Navier Stokes equations. Table 10 displays the results obtained.
We see that access data essential, as training solely on the PDE residual leads to high data loss. Also, training with data on fixed grid results in poor generalization to higher resolutions. By contrast, with Autograd, small number of randomly sampled points is sufficient to get excellent results. This demonstrates that random sampling allows the Autograd mGNO FNO to capture correct underlying physics instoad of overfitting to a fixed grid. Also, including a physics loss is critical. Models trained only on data loss exhibit extremely large PDF residuals, while incorporating the physics term, even at low resolution with randomly sampled points, ensures good PDE fidelity. This indicates that the physics loss guides the network toward physically consistent solutions, and very high-resolution PDE evaluation is not always necessary.
Finally, none of the results obtained with finite differences satisfactory. It worth noting that these trained with physics loss fixed grid (given the added complexity and reduced reliability of non-uniform finite differences), which can limit generalization across resolutions. However, even accounting for this, using finite differences physics losses led to results significantly worse than training purely on data, both in terms of PDE residual and data loss at the test resolutions. This reveals a clear mismatch between the finite-difforence PDE loss and the reference data, suggesting that the finite-iifference derivatives were not sufficiently accurate the training resolutions, even when using as many as 50,000 randomly sampled points.
19

---

Published in Transactions on Machine Learning Research (10/2025)
Table 10: Results of training mGNO FNO models with hybrid loss for the Navier Stokes equations, using varying resolutions for the physics term (PDE Res.) and the data term (Data Res.) during training. When the resolution is denoted by single number N, it indicates that N points are randomly sampled at each training iteration. In contrast, when the resolution is written as N₂ x N₂ x Nr, it refers to using a fixed grid with the specified spatial and temporal resolutions. The table reports both PDE residuals and data losses at the training resolutions, and at the testing resolutions (fixed to 32 x 32 16 for the PDE residual and 64 64 50 for the data loss). Derivatives are computed with Autograd and finite differences. Cells highlighted indicate suboptimal results (orange) or poor generalization (red) at the test resolution.
AUTOGRAD
PDE Res.
Data Res.
PDE (Train)
PDE (Test)
Data (Train)
Data (Test)
10000
32
50
3.25 10-3
6.99 10-5
2.91 10-4
1.50 10-2
10000
50000
1071.10-9
2.13 10
1.06 10
1.05 10
10000
10000
1052-10-
2.04-10-5
1.11
1.09 10
10000
1000
2
3.11-10-5
1.53-10-
2.27 10-4
10000
100
1.70-10-4
1
1.32 10-3
1.87 10
10000
10-
10-7
2.04
1000
32
32 50
2074-10-
6.28-10-7
1.57 1 10-5
1.-0-1--
1000
50000
2.48.10
10-7
5.51 1.511105
5.79. 10
1000
10000
30-
60-2-10-7
5.44
5.28 10
1000
1000
6.74-10-6
2.02-10-1
1.48 10-
1.31 13-
1000
100
-
3.55 10-3
3.95 10-3
1000
4
1.02
2.84 10
100
32
32 50
-
6.58 10:8-00- 10-6
-
100
50000
10
5
3.49 3
3.53 10 10
100
10000
1
2.47-10-6
-
8.06
100
1000
1.97 10-6
2074-10-6
8.59 10-5
7.69 10-5
100
100
6
-
2.50 10-4
2.32 10-4
100
2.92.
-
2.01 2.01 10
25
32
32 x 50
9.28 10-3
1.83 10-7
8.04 10-6
1.50 10
25
50000
9.98
9
6.63 10-5
6.84 10
25
10000
20
1.07 -
1.05
25
1000
1.27 1.27
1.63 1663-10-
1.16 10 10-4
1.35 105-100 10-4
25
100
2.46 20-5
3.53 10-5
2
3.08 10-4
25
1.96 10-7
1.02 10-6
1.92 10
FINITE DIFFERENCES
PDE Res.
Data Res.
PDE (Train)
PDE (Test)
Data (Train)
Data (Test)
32
32
25
32 x 32 50
9.14 10
1.09 106
3.07 x 10
1.29 10
16
16
25
32 x 32 50
5.16 10-4
3.42 105
1.78 X 10-3
4.28 x 10-2
16
x
16
x
12
32 x 32 x 50
1.76 x 10
2.45 x 107
8.28 x 10
1.12 x 10
32
32
25
50000
3.58 10
1.93 x 105
3.43 x 10-3
3.51 X 10
32
32
25
10000
3.20 10-3
8.95 104
3.73
-0-3
3.54 -0-3
NO PHYSICS LOSS
PDE Res.
Data Res.
PDE (Train)
PDE (Test)
Data (Train)
Data (Test)
32 x 32 50
3.55 10a
5.30 10
1.56 10
50000
2.26 104
1.87 10
1.82 10
10000
2.83 104
1.97 10 10-6
1.91 10 10-6
1000
2.47 105
3.84 10
3.51 10
100
2.83 104
4.18 10
6.54 -
20

---

Published in Transactions on Machine Learning Research (10/2025)
4.5
Airfoil Inverse Design
4.5.1
Description of the Forward Mapping
We consider the transonic flow over an airfoil (ignoring the viscous elfect), governed by the Euler equation,
Op
+ v (p'v) 0,
(22)
dt
OF
V (EE + +pv) = 0,
(23)
0p'v v
+ (p 2 v + pl) O,
(24)
where of is the fluid density, is the velocity vector, p is the pressure, and E is the total energy. The far-field boundary condition is P:o Pxe = 1, Mx 0.8, АoA 0, where Ms0 is the Mach number and AoA is the angle of attack, and the no-penetration condition is imposed at the airfoil.
We use the same dataset as Li et al. (2022). where the shape parameterization of the airfoil follows the design element approach (Farin| 2014). The initial NACA-0012 shape is mapped onto a 'cubic' design element with control nodes in the vertical direction with prior d U[-0.05, 0.05). That initial shape morphed to a different shape following the displacement field of the control nodes. The dataset contains 1000 training samples and 200 test samples generated using second-order implicit finite volume solver. The C-grid mesh with (220 x 50) quadrilateral elements is used and adapted near the airfoil, but not around the shock.
For the forward pass, the mesh point locations and signed distance functions are passed through a trained differentiable mollified GINO model of the form
GmGINO
Geecoder
GeNO
Geneoder
(25)
to produce an approximation of the pressure field p.
4.5.2 Inverse Design
For the inverse problem, parametrize the shape of the airfoil by the vertical displacements of few spline nodes, and set the design goal to minimize the drag-lift ratio. The parametrized displacements of the spline nodes are mapped to a mesh which is passed through the mGINO to obtain a pressure field, from which we can obtain the drag-lift ratio. We optimize the vertical displacement spline nodes by differentiating through this entire procedure. depiction of the airfoil design problem is given Figure 13
Backpropagate gradients to update control
mGINO
drag ratio
cantral spline nodes
point cloud
pressure p
The airfoil design problem. Parametrized vertical displacements control spline nodes generate point cloud, which passed through the differentiable mGNNO to obtain a pressure field p, from which we can compute the diag-lift ratio. We update the control nodes minimize the drag-lift ratio by differentiating through this entire procedure.
Figure 13:
21

---

Published in Transactions on Machine Learning Research (10/2025)
As result of this optimization process, obtain an airfoil design with drag coefficient of 0.0216 and lift coefficient of 0.2371, based on the model prediction. Using the numerical solver on this optimal design, we verify these predictions and obtain a similar drag coefficient of 0.0217 and a similar lift coefficient of 0.2532.
This yields a drag-lift ratio of around 0.09, outperforming the optimal drag-lift ratio of 0.14 reported by Li et al. )022) (drag 0.04 and lift 0.29, obtained using a Geo-FNO instead of mGINO).
Discussion
We proposed mGNO, fully differentiable version of GNO, allow for the use of automatic differentiation when computing derivatives, and embedded within GINO to learn efficiently solution operators of families of large- scale PDEs with varying geometries without data. The proposed approach circumvents the computational limitations of traditional solvers, the heavy data requirement of fully data-driven approaches, and the generalization issues of PINNs.
The use of a physics loss proved critical and was sufficient to achieve very good results with data at a very low resolution and even in the absence of data, demonstrating the data efficiency of the physics-informed paradigm. This highlights the need for elficient and accurate methods to compute physics losses, improve data efficiency, and regularize neural operators. Autograd can compute exact derivatives in a single pass seamlessly across complex geometries and enables higher-order derivatives with minimal additional cost. However, it can be memory-intensive for deep models, possibly making it less efficient than finite differences for simpler problems. Despile these limitations, its accuracy, capability to handle complex geometries, and scalability for complicated learning tasks make it the preferred differentiation method in PINNs, and a promising approach in our plysics-informed neural operator setting on complex domains.
Using Autograd instead of finite differences led to 20x reduction of the relative L2 data loss for Burgers' equation regular grids, suggesting that the Autograd physics loss better captured the physics underlying the data. In hybrid training with noisy reference data, mGNO o FNO remained robust and accurate, maintaining low PDE residuals and data loss. It also excelled at learning the lid cavity flow example, where a small number of randomly sampled points proved sufficient with Autograd to achieve excellent. results. Here again, including a physics loss, potentially at a small number of randomly sampled points, proved crilical. Autograd mGINOO performed seamlessly for the Poisson and hyperelasticity equations on unstructured point clouds, while finite differences were not sufficienlly accurate at the training resolution used and would need at lcast 9x more points to compute reasonable derivatives. Autograd mGINO achieved relative error 2-3 orders of magnilude lower than the Meta-PDE baselines for comparable running time, and enjoyed speedups of 20-25x and 3000-4000x compared to the solver for similar accuracy on the Poisson and hyperelasticity equations. We also demonstrated with an airfoil design problem that mGINO can be used scamlessly for solving inverse design and shape optimization problems on complex geometries.
Future Directions
Our framework lays the groundwork for further transformative advances in physics-informed machine learning.
Integrating techniques that have proven effective in PINNs and PINOs, such as adaptive loss reweighting and targeted sampling, could dramatically enhance training efficiency and solution accuracy, while advanced
automatic differentiation strategies like the Stochastic Taylor Derivative Estimator (Shi et al. 2024 can help make high-dimensional PDEs with exact derivatives more computationally tractable. Learning data-driven mollifiers directly from PDE residuals could allow the model to adaptively tailor its smoothing properties to the structure of the underlying PDE, potentially improving both stability and generalization across regimes. Adaptive meshing further leverages the model's flexibility with arbitrary point clouds, concentrating computational effort on regions with steep gradients, singularities, or other complex dynamics, thereby accelerating training and enabling highly efficient forward and inverse design.
22

---

Published in Transactions on Machine Learning Research (10/2025)
Data and training efficiency can be further improved by leveraging multiscale strategies. As shown throughout the paper, mGNO achieves high accuracy and strong data efficiency by combining low-resolution data with high-resolution physics losses. This aligus with growing recent evidence for the benefits of multiscale methods. For example, Ahamed et al. (2025) introduces a CNN training scheme across spatial resolutions, where gradual refinement. maintains stable gradients and improves efficiency and generalization when direct. high-resolution training would be too costly. Gal ct al. (2025) propose a coarsc-to-fine GNN framework using graph coarsening and multiscale gradients to reduce memory and computation without sacrificing accuracy. While we do not apply hierarchical coarsening, subsampling PDE residual points similarly reduces cost with minimal accuracy loss and suggests potential for adaptive sampling. George et al. (2024) extends the idea to FNOs with progressive training in spatial resolution and spectral modes, lowering cost, improving generalization, and producing more compact models than full-resolution training. This highlights the broader advantage resolution adaptivity, which mGINO can exploit when handling point clouds, for example, by gradually incrcasing Fourier modes alongside resolution during training. Together, these results underscore the value of exploring multiscale extensions of mGNO, such as mollified kernels of varying radii, enabling simultaneous capture of short- and range interactions and addressing multiscale and nonlocal phenomena that challenge conventional architectures.
Together, these directions underscore the transformative potential of our approach for modeling complex, high-dimensional PDFs, opening pathways for both predictive accuracy and computational efficiency that extend well beyond current state-of-the-at methods.
Acknowledgments
Ryan Lin supported by the Caltech Summer Undergraduate Research Fellowships (SURF) program. Anima Anandkumar supported in part by the Bren endowed chair, ONR (MURI grant N00014-23-1-2654), and the A12050 Senior Fellow Program at Schmidt Sciences.
References
Shadab Ahamed, Niloufar Zakariaci, Eldad Haber, and Moshe Eliasof. Multiscale training of convolutional neural networks, 2025. URL
org/abs/25011...... 12739
Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie Rognes, and Garth Wells. The fenics project version 1.5. 3, 01 2015. doi: 10.11588/ans.211511000. 100.0553..
Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi Li, Miguel Liu-Schiaffni, Jean Kossaifi, and Anima Anandkumar. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, pp. 1 9, 2024.
Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimencz Rezende, et al. Interaction networks for learning about objects, relations and physics. Advances in neural information processing systems, 29, 2016.
Atilm Günes Baydin, Barak A. Pearlmutter, Alexcy Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: a survey, J. Mach. Learn. Res., 18(1):5595 5637, January 2017. ISSN
1532-4435.
Julius Berner, Miguel Liu-Schiaffini, Jean Kossaifi, Valentin Duruisseaux, Boris Bonev, Kamyar Azizzade- nesheli, and Anima Anandkumar. Principled approaches for extending neural architectures to function spaces for operator learning. 2025. URL httpp::/aaxxv..or//aaaee
10973.
Rafael Bischof and Michael Kraus. Multi-objective loss balancing for physics-informed deep learning. 2021.
J. W. Burby, Q. Tang, and R. Manlik. Fast neural Poincaré maps for toroidal magnetic fields. Plasma Physics and Controlled Fusion, 63(2)0024001, dec 2020.
23

---

Published in Transactions on Machine Learning Research (10/2025)
Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George E. Karniadakis. Physics- infermedd neural networks (PINNs) for fluid mechanics: review. Acta Mechanico Sinica, 37, 01 2022. doi: 101100774110088022--144441
Nithin Chalapathi, Yiheng Du, and Aditi S. Krishnapriyan. Scaling physics-informed hard constraints with mixture-of-cxperts. In The Twelfth International Conference on Learning Representations, 2024.
Renyi Chen and Molei Tao. Data-driven prediction of general hamiltonian dynamics via learning exactly- symplectic maps. In Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 1717-1727. PMLR, 18-24 Jul 2021.
Zhao Chen, Vijay Badrinarayanan, Chen Yu Lee, and Andrew Rabinovich. GradNorm: Gradient normalization for adaptive loss balancing in deep multitask networks. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedlings of Machine Learning Research, pp- 794-803. PMLR, 2018.
Valentin Duruisseaux, Joshua W. Burby, and Qi Tang. Approxination of nearly-periodic symplectic maps via structure preservvng neural networks. Scientific Reports, 13(8351), 2023.
Valentin Duruisscaux, Miguel Liu-Schiaffni, Julius Berner, and Anima Anandkumar. Towards enforcing hard physics constraints in operator learning frameworks. In TCMT, 2024 AT for Science Workshop, 2024.
Lawrence C. Evans. Partial differential equations. American Mathematical Society, Providence, R.I., 2010.
Gerald Farin. Curves and surfaces for computer-aided geometric design: a practical guide. Elsevier, 2014.
Chelsea Finn, P. Abbeel, and Sergey Levine. Model-agnostic meta-learning for fast adaptation of deep networks. In International Conference on Machine Learning, 2017.
Scbastian Flennerhag, Pablo G. Moreno, Neil D. Lawrence, and Andreas Damianou. Transferring knowledge
across learning processes, 2019.
Eshed Gal, Moshe Eliasof, Carola-Bibiana Schönlieb, Ivan 1. Kyrchei, Eldad Haber, and Eran Treister.
Towards efficient training of graph neural networks: A multiscale approach, 2025. URL htttps://axii>
org/abs/2503. 19666
Adarsh Ganeshram, Haydn Maust, Valentin Duruisseaux, Zongyi Li, Yixuan Wang, Daniel Leibovici, Oscar Bruno, Thomas Hou, and Anima Anandkumar. Fe-pino: High precision physics- informed neural operators via fourier contimuation, 2025. URL jtttps://arxv. a2211 15960
Robert Joseph George, Jiawei Zhao, Jean Kossaifi, Zongyi Li, and Anima Anandkumar. Incremental spatial and spectral learning of neural operators for solving largo-scale PDEs Transactions on Machine Learning Research, 2024. ISSN 2835-8856. URL httts:://pperreviet nott/formm7dxxx1ePP00pp
Luca Ghafourpour, Valentin Duruisseaux, Bahareh Tolooshams, Philip H. Wong, Costas A. Anastassiou, and Anima Anandkumar. Noble noural operator with biologically latent emboddings to capture experimental variability in biological neuron models, 2025. URL https: www.ziiv.rr//a//2555552 04536
Vignesh Gopakumar, Stanislas Pamela, Lorenzo Zanisi, Zongyi Li, Anima Anandkumar, and MAST Team. Fourier neural operator for plasma modelling. arXiv preprint arXiv:2302.00 06542, 2023.
Paula Harder, Alex Garcan, Venkatesh Ramesh, Qidong Yang. Prasanna Sattigeri, Daniela Szwareman, Campbell Watson, and David Rolnick. Hard-constrained deep learning for climate downscaling, 2024.
A. Ali Heydari, Craig Thompson, and Asif Mehmood. SoftAdapt: Techniques for adaptive loss weighting of neural networks with multi-part loss functions. 2019.
Ameya D. Jagtap, Ehsan Kharazmi, and George E. Karniadakis. Conservative informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. Computer Methods in Applied Mechanics and Engineering, 365:113028, 2020. ISSN 0045-7825. doi: 00.1016/..ma..200.111 13028.
24

---

Published in Transactions on Machine Learning Research (10/2025)
Chiyu "Max" Jiang, Karthik Kashinath, Prabhat, and Philip Marcus. Enforcing physical constraints in CNNs through differentiable PDE layer. In ICLR 2020 Workshop on Integration of Deep Newral Models and Differential Equations, 2020.
P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis. SympNets:
Intrinsic peructure-preserving
symplectic networks for identifying Hamiltonian systems. Neural Networks, 132(C), 12 2020.
G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Plysics-informed machine learning. Nature Reviews Physics, 3(6):422 440, 2021.
Jean Kossaifi, Nikola Kovachki, Zongyi Li, David Pitt, Miguel Liu-Schiaffnn, Valentin Duruisseaux, Robert Joseph George, Boris Bonev, Kamyar Azizzadenesheli, Julius Berner, and Anima Anandkumar. A library for learning neural operators, 2025. URL https://aaxxt wwtxxiiv.rrg//oooo abs abs//aa/....////// 2412. 10354
Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for Fourier neural operators. J. Mach. Learn. Res., 22(1), 2021. ISSN 1532-4435.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1-9 97, 2023.
Aditi S. Krishnapriyan, Amir Gholami, Shandian Zhe, Robert M. Kirby, and Michael W. Mahoney. Char- acterizing possible failure modes physics-informed neural networks. In Neural Information Processing Systems, 2021.
Thorsten Kurth, Shashank Subramanian, Peter Harrington, Jaidoep Pathak, Morteza Mardani, David Hall, Andrea Miele, Karthik Kashinath, and Animashree Anandkumar. FourCast Accelerating global high-resolution weather forecasting using adaptive Fourier neural operators. 2022. doi: 10.48550/ark arvv.. 2208.05419.
Raphael Leiteritz and Dirk Pflüger. How to avoid trivial solutions in physics-informed neural networks. ArXiv, abs/ (bls/2122.0560, 2021.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint ar Xiv:20010. 08895, 2020a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kanshik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Graph kernel network for partial differential equations. arXie preprint arXi::2003.034 03345, 2020b.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkunar, et al. Fourier neural operator for parametric partial differential equations. In ICLR, 2021a.
Zongyi Ti, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Axizzadenesheli, and Anima Anandkumar. ACM/JMS Journal of Data Science, 2021b.
neural operator for learning partial differential equations.
Zongyi Li, Danicl Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with
learned deformations for pdes general geometries. arXie preprint arXiv:2207. 052009, 2022.
Zongyi Li, Nikola Borislavov Kovachki, Chris Choy, Boyi Li, Jean Kossaifi, Shourya Prakash Otta, Moham- mad Amin Nabian, Maximilian Stadler, Christian Hundt, Kamyar Azizzadenesheli, and Anima Anandkumar. Geometry-informeed neural operator for large-scale 3d pdes. arXiv preprint X3W:2339.00583 2023.
G.R. Liu. Meshfree Methods: Moving Beyond the Finite Element Method. CRC Press, second edition edition, 2009. doi: 10.1 12200/97714220 42000220081
25

---

Published in Transactions on Machine Learning Research (10/2025)
Anders Logg, Garth Wells, and Kent-Andre Mardal. Automated Solution of Differential Equations by the Finile Flement Method: The FFniCS Book, volume 84. 04 2011. ISBN 978-3-642-23098-11 doi: 10.1007/978 0717779978---22-280 23099-8.
Huakun Luo, Haixu Wu, Hang Zhou, Lanxiang Xing, Yichen Di, Jianmin Wang, and Mingsheng Long. +: An accurate neural solver for pdes on million-scale geometries, 2025. URL https://arxix: org/ 2tg/abs/2502..0445
Dougal Maclaurin, David Duvenaud, and Ryan P Adams. Autograd: Effortless gradients in numpy. In ICML 2015 AutoML Workshop, 2015.
Qinglu Meng, Yingguang Li, Xu Liu, Gengxiang Chen, and Xiaozhong Hao. A novel physics-informed neural operator for thermochemical curing analysis of rrrbonfibree-riiforredd thermosetting composites. Composite Structures, pp. 117197, 2023.
Arvind T. Mohan, Nicholas Lubbers, Misha Chertkov, and Daniel Livescu. Embedding hard physical constraints neural network coarso-graining of threc-dimensional turbulence. Phys. Rev. Fluids, 8:014604, Jan 2023. doi: 10.1103/ 10.1103/PhysRev Physseeeluidss..00 4.04
Geoffrey Negiar, Michael W. Mahoncy, and Aditi S. Krishnapriyan. Learning differentiable solvers for systems with hard constraints. 2022.
Johannes T.B. Overvelde and Katia Bertoldi. Relating pore shape to the non-linear response of periodic elastomeric structures. Journal of the Mechanics and Physics of Solids, 64:351-300 336, 2014. ISSN 0022-5096. doi: https://doi.orr, 10. 10.1000.1...... htw.:/...pp..201...1...
Tian Qin, Alex Beatson, Deniz Oktay, Nick McGreivy, and Ryan P. Adams. Meta-PDF: Learning to Solve PDEs Quickly Without a Mesh. 2022.
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 Compulational Physics,
707, 2019. ISSN 0021-9991.
doi: 10. 10.1016/1jee:2010810000
Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. ArXiv, abs/1711.10561, 2017a.
Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations. ArXie, abs/1711.10566, 2017b.
Jack Richter-Powell, Yaron Lipman, and Ricky T. Q. Chen. Neural conservation laws: A divergence-free perspective. In Advances in Neural Information Processing Systems, 2022.
Shawn G Rosolsky, Hani AI Majed, and FA Huerta. Applications of physics informed neural operators. Machine Learning: Science and Technology, 4(2):025022, 2023.
Zekun Shi, Zheyuan Hu, Min Lin, and Kenji Kawaguchi. Stochastic taylor derivative estimator: Efficientt amortization for arbitrary differential operators. In The Thirty-eighth Annual Conference Newral Information Processing Systems, 2024.
Yuchen Song, Danshi Wang, Qirui Fan, Xiaotian Jiang, Xiao Luo, and Min Zhang, Physics-informed neural operator for fast and scalable optical fiber channel modelling in multi-span transmission. In 2022 European Conference on Optical Communication (ECOC), pp. 1-4. IEEE, 2022.
Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. SIAM Journal on Scientific Computing, 43(5):33055 A3081, 2021. doi: 10.1137/20M1311044.
Gege Wen, Zongyi Li, Qirui Long, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M. Benson. Real-time high-resolution CO2 geological storage prediction using nested Fourier neural operators. Energy Enairon. Sci., 1741, 2023. doi: 10. 10.1039 F0039/D2FF04204F.
26

---

Published in Transactions on Machine Learning Research (10/2025)
Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long. Transolver: A fast transformer solver for pdes on general geometries, 2024. URL https://arxxv..rrrt abaaaaaaaaa 22022 02366
Lanxiang Xing, Haixu Wu, Yuezhou Ma, Jianmin Wang, and Mingsheng Long. Helmfluid: Learning helmholtz dynamics for interpretable fluid prediction. In International Conference Machine Learning, 2024.
Zhanhong Ye, Xiang Huang, Leheng Chen, Hongsheng Liu, Zidong Wang, and Bin Dong. Pdeformer: Towards a foundation model for one-dimensional partial differential equations, 2025. URL https://arxiv.orr abs/2402. 12652
Jeremy Yu, Lu Lu, Xuhui Meng, and George E. Karniadakis. Gradient-enhanced physics-informee neural nelworks for forward and inverse PDE problems. Compuler Methods in Applied Mechanics and Engineering, 393: 19422,,, 2022. ISSN 0045-7825. doi: 10.1016/j.cma. 20...m..22........
Hang Zhou, Yuezhou Ma, Haixu Wu, Haowen Wang, and Mingsheng Long. Unisolver: Pde-conditional transformers are universal pde solvers, 2025. URL http::///aaxii..orr//////22000
27

---

Published in Transactions on Machine Learning Research (10/2025)
A
FNO and GINO Architecture Diagrams
We first depict the Fourier Neural Operator (FNO)
[Li et al.
2020a architecture:
a(x)
P
Fourier layer
Fourier layer 2
Fourier layer T
Q
u(x)
Fourier layer
F
R
v(x)
W
Figure 14: The Fourier Neural Operator (FNO) architecture (extracted from Li et al. (2020a
Next, depict the Geomelry-Informee Neural Operator (GINO)
Li el al. 2023 architecture:
Input Samples
discretize
GNO
FNO
GNO
query
Prediction
Operates on input geometry
Operates on (latent) regular geometry
Projects to output geometry
MAAAA WMMN
FNO layers
Figure 15: The Geometry- Informmd Neural Operator (GINO) architecture (extracted from Li et al. 2023
B
Code
The PyTorch code used for the experiments presented in this paper has been added to the open-source neural operator library from Kossaifi et al. (025) at http://hhh...omooeeee////e///eer
28

---

Published in Transactions on Machine Learning Research (10/2025)
C Choices of Hyperparameters
For all experiments, the "optimal' hyperparameters used (including weight functions and loss coefficients) were obtained by conducting grid search on subset of hyperparameters.
For the radius cutoff used in the weight functions of mGNOs, we use the same radius as for the GNO to avoid the additional cost of carrying a second neighbor search.
We use the following Mean Squared Error (MSE), Relative Squared Error, L2 error, and Relative L2 errors metries in experiments:
MSF (ypreds Ytrue)
Σ
(6) pred
(1) Virue
(t] Virue
1 pped
RELATIVESSQARREEEERRROKO.... ytrue)
Btrue
E
N
L.2 (.2 Ptrue)
C
Σ
(Upredt
Virue
CVEM" Upreet
C.
(6) Berue
( Ptrue)
CVER
Virue
+ 8
where is small positive number for numerical stability, and the constant C is scaling constant taking into account the mcasure and dimensions of the data, to ensure that the loss is averaged correctly across the spatial dimensions.
C.1 Burgers' Equation
For Burger's equation, the input initial condition un(z) given on a regular spatial grid in the domain [0, 1]2 is first duplicated along the temporal dimension obtain 2D regular grid of resolution 128 x 26, and then passed through a 2D FNO and mollified GNO to produce predicted function
(GmGNO Grno)(u))
approximating the solution u(x, l).
For this experiment, the 2D FNO has 4 layers, each with 26 hidden channels and (24, 24) Fourier modes, and used Tucker factorization with rank 0.6 of the weights. The mGNO uses the half cos weight function with radius of 0.1, and 2-layer MLP with [64, 64) nodes.
The resulting model has 1, 019, 569 trainable parameters, and was trained in PyTorrh for 10, 000 epochs using the Adam optimizer with lcarning rate 0.002 and weight docay 10-6, and the ReduceLRRnPlatcan scheduler with factor 0.9 and patience 50.
29

---

Published in Transactions on Machine Learning Research (10/2025)
C.2 Nonlinear Poisson Equation
For the nonlinear Poisson equation, the mesh coordinates within [-1.4, 1.4]², signed distance functions, source terms, and boundary conditions, are passed through
GMIGNO GeNO 9GNO
model to produce an approximation to the solution u.
For this experiment, the imput GNO has a radius of 0.16, and 3-layer MLP with [256, 512, 256] nodes. The 2D FNO has 4 layers, cach with 64 hidden channels and (20, 20) Fourier modes, and acts on a latent space of resolution 64 x 64. The output mGNO uses the half cos weight function with a radius of 0.175, and a 3-layer MLP with [512, 1024, 512] nodes.
The resulting model has 8, 691, 972 trainable parameters, and was trained in PyTorch for 300 epochs using the Adam optimizer with learning rate 0.0001 and weight decay 10-6, and the ReduceLLROnllaten scheduler with factor 0.9 and patience 2. We used 7000 samples for training and 3000 samples for testing.
C.3 Navier-Stokes Equations
For the Navier Stokes equations, employ a physics-informed neural operator approach that learns the mapping from spatial-temporal coordinates (x, v. to the velocity-presse pressure field eu,,,)). The model architecture combines a 3D FNO with a mollified GNO to produce a predicted function
= (ImGNo GFNo) v.
approximating the solution (u(2, y, t), f(x, y, t) p(z,))p((2,,) p(z, t)).
The 3D FNO has 3 layers, each with 16 hidden channels and (16, 16, 16) Fourier modes, and uses Tucker factorization with rank 0.1 of the weights. The mGNO uses the half cos weight function with radius 0.016, a 2-layer channel MLP with [128, 64) nodes. The resulting model has 431, 183 trainable parameters, and was trained in PyTorc for 15, 001 epochs using the Adam optimizer with learning rate 0.001 and the RedueeLLROnPlateal scheduler with factor 0.8 and patience 200.
The training loss combines (up to) four loss components: (1) PDE residual enforcing the Navier-Stoker equations varying numbers of randomly subsampled points from a high-resolution 256 x 256 80 grid, (2) data loss matching velocity predictions to target data on either varying number of points subsampled from a 64 x 64 x 50 grid or a static 32 x 32 x 50 grid, (3) initial condition loss enforcing zero initial velocity, and (4) boundary condition loss enforcing no-slip conditions on walls and lid-driven flow on the top boundary. Loss weighting was handled automatically and adaptively using ReLoBRaLo (Bischof & Kraus 2021).
C.4 Hyperelasticity Equation
For the hyperelasticity equation, the mesh coordinates within (0, 1/2 and signed distance functions are passed through
Glecoder o GeNO gencoder
model to produce an approximation to the solution u the hyperelasticity equations.
For this experiment, the input GNO has radius of 0.05625, and 3-layer MLP with [128, 256, 128] nodes. The 2D FNO has 4 layers, cach with 64 hidden channels and (20, 20) Fourier modes, and acts on a latent space of resolution 32 x 32, The output mGNO uses the half cos weight function with a radius of 0.1125, and a 3-layer MLP with [1024, 2048, 1024) nodes.
The resulting model has 11, 67,, 211 trainable parameters, and was trained in PyTorc for 400 epochs using the Adam optimizer with learning rate 0.0001 and weight decay 10-6, and the PedneeLROnPlateau scheduler with factor 0.9 and patience 10. We used 1000 samples for training and 1000 samples for testing.
30

---

Published in Transactions on Machine Learning Research (10/2025)
C.5 Airfoil Inverse Design
For the airfoil design forward problem, the 220 x 50 mesh point locations within [-40, 40j² and signed distance functions are passed through differentiable mollified GINO model of the form
IMGINO = mmNNO GrNO MGGNO
produce an approximation of the pressure field p.
For this experiment, the input mGNO uses the half weight. function with radius of 2, and 3-layer MLP with [128, 128, 128] nodes. The 2D FNO has 3 layers, each with 16 hidden channels and (36, 36) Fourier modes, and acts a latent space of resolution 64 x 64. The output mGNO uses the half colf cos weight function
with a radius of 6, and a 3-layer MLP with [128, 128, 128] nodes.
The resulting model has 1, 162, 546 trainable parameters, and was trained in PyTorch for 750 epochs using the Adam optimizer with learning rate 0.0001 and weight docay 10-9, 09 and the RodueeLLROnPlatea scheduler with factor 9 and patience 5. We used 1000 samples for training and 200 samples for testing.
D GNOO Layer Pseudocode
Here, we provide a simplified example of Py Torch pseudocode for the mGNO layer
GmcNo(U)(x)= = J/-
v)v(y) dy,
with the half cos weight function Whalf cos (I, y) = (z) 005 0.5))s(=d) where d = +| - for--112/v² /112/12.
def half_com cos weight Auncioondistt, radius=1., scale=1.);
_oiighh
return scale (0.5 torch. cos (torch. pi diata**2 / radius**2) 0.
def mGNO_layer
V: Tensor [bs, n_in, codim], # discretization of the function to transform coordinates of the discretization of query locations
y: Tensor [bs, m_in, dim]
x: Tensor [bs, n_out, dim],
delta: Tensor [bs, n_in], radius=None,
quadrature weights when approzimating integral # radius of the mollified GNO
10
weight_fn=None
weight function w,
net: torch. mn Module # neural network parometrizing the kernel, e.g. MLP Tensor [bs, n_out, codim]:
10
14
Kernel evaluation shape = [bs, n_in, n_out, dim]
10
kernel_ing [y. unsqueeze (2) expand (shape), x. unsqueeze (1) expand (shape)) kernel net (torch.. cat (kernel_ing dim=-1))
18
# Weighted aggregation using the quadrature weights output = kernel v. unsqueeze (1) delta. delta. vbaw m_in, 1, 1)
20
# Hollefication radius is not None:
23
dists = cdist(y, x) get distances between input and query locations
25
output [dists radius, :] = give weight outside ball of radius r if veight_fa fn is not None:
27
Apply the weight function
output output weight_ fn (dists, radius). unsqueeze (-1)
20
return output, sum (dim=-2)
31

---

Published in Transactions on Machine Learning Research (10/2025)
E Finite Differences on Point Clouds
On a regular grid, standard well-known stencil formulas are available to compute first-order and higher-order derivatives using finite differences. However, on arbitrary point clouds, the varying distances between points have to be taken into account, and a different stencil needed for each point at which the derivatives need to be computed. We detail below one strategy to obtain these stencil formulas in 2D.
We consider the case where have arbitrary point 2D cloud with N points, {(2,, y.) and suppose that the function values (f(x; 8s)} of function f are known at these points. The goal is to approximate partial derivatives of f at any other point (2,9) in the domain. We start with first-order derivatives, and write them as a finite difference with unknown stencil coefficients:
0f dx ( 1)
(x) (a,,)i))
0f dy
(t. 1)
(v) ((),
(26)
1=1
To find the stencil coefficients c and o enforce that the approximation holds exactly for the functions 1, and that is, we enforce that the approximation holds true for any polynomial of degree 1 in 2D. This results in the following systems of equations for the stencil coefficients e( and c
(2)
(v)
Σ
Σ
(a: 8) 1,
Σ
(0) (v
(27)
and
(v) = 0.
(2) (#: F) =0
Σ
(v) (v: - 5) =
(28)
t=1
i=1
t=1
This can be written as
xx =b(²),
Ac b(u),
(29)
where
1
A
T1 T.
D22 5
EN $
(30)
91
1/2
UN
and
c(x)_
(1))
b
b =
(31)
These systems of equations can be solved using least squares:
(z) (ATA)-'ATbb
((u) =
(32)
Plugging these coefficients in Equation
26 gives the desired approximation to the first-order derivatives at f.
We emphasize that this procedure needs to be repeated for every point (ž, Ø) at which the derivatives need to be evaluated, since the location of the point affects the entries of the matrix A. Note that it is not necessary and not recommended to use all points to approximate the derivatives, and one should instead identify a subset of nearest neighbors from which the finite differences can be computed.
To compute higher-order derivatives, one could follow similar strategy (which would lead to complicated system of equations to solve), first evaluate first derivatives on the point cloud and repeat the above procedure by replacing of by its appropriate partial derivative.
32

---

Published in Transactions on Machine Learning Research (10/2025)
F
Training Times of 2D mGNOs with Autograd and Finite Differences
We compare the training times per epoch when training 2D mGNOs with Proisson using Autograd, finite differences on uniform grids, and finite differences on non-uniform (NU) grids. Figures 16 and 17 display the training times per epoch versus the latent resolution and output resolution, respectively.
As expected, we see from Figure 16 that Autograd is more expensive than finite differences on a uniform grid fixed latent and output resolutions. We can also see, by looking at individual columns corresponding to fixed latent space resolutions (i.e. the same model architecture), how finite differences at higher output resolutions compare to Autograd at lower output resolutions. In particular, non-uniform finite differences are often more expensive than Autograd at a slightly lower resolution. Figure 17 quantifies how running time increases when the latent space resolution of the model increases.
Finite Differences
NU Finite Differences
Autograd
128
64
160
160
128
0.5
64
128
96
128
128
128
96
500
256
128
Time per epoch (s)
0.2
64
96
96
64
32
64
488
0.1
128
48
64
48
32
256
32
128
64
48
48
32
64
32
16
8
0.05
16
16
32
16
16
16
32
32
64
8x8
16x16
32x32
48x48
Latent Resolution
Figure 16: Time per epoch when training 2D mGNOs with Lroison using Autograd, finite differences on uniform grids, and finite differences on non-uniform (NU) grids. We display the training times per epoch different latent space resolutions (x-oris) and different output resolutions (the numbers next to the points denote the number of points along each dimension).
33

---

Published in Transactions on Machine Learning Research (10/2025)
Finite Differences
32
64
16
16
32
16
NU Finite Differences
0.4
Autograd
Time per epoch (s)
0.2
0.1
0.02 8x8
16x16
32x32
64x64 Output Query Resolution
128x128
256x256
Figure 17: Time per epoch when training 2D mGNOs with Poosson using Aulograd, finite differences on uniform grids, and finite differences on non-umiform (NU) grids. We display the training times per epoch at different output resolutions (x-oris) and different latent space resolutions (the numbers at the end of the curves denote the number of points along each dimension).
G Automatic Differentiation with Non-Differentiable Components
In our experiments, the entire mGNO and mGINO architectures are designed such that all of their components are differentiable, ensuring that derivatives of any order can computed correctly and accurately using automatic differentiation. For instance, we avoid non-differentiable activation functions such as ReLU (Rectified Linear Unit), which introduce discontinuities in their derívatives, and instead smooth alternatives like GeLU (Gaussian Error Lincar Unit). All the other operations in the network, including the Fourier layers, pointwise linear neural networks, and pointwise linear operators, are also implemented in differentiable manner to avoid issues related to non-differentiabiliy:
While automatic differentiation is naturally suited to smooth functions, can accommodate non-diffcrentiable components. For example, ReLU non-differentiable at 0, yel. Autograd compules exact. gradients on each side and assigns a subgradient at 0 (e.g. or 1). This works well in practice since exact zero inputs are rare and subgradient- optimization remains effective. Other strategics exist for non-smooth components. Clarke subderivatives formalize generalized derivatives non-smooth points. Proximal methods split updates into smooth gradient step and proximal mapping for the non-smooth part. Randomized smoothing adds small perturbations to inputs and averages outputs. Straight Through Estimators use the true non-differentiab function the forward pass while replacing its gradient in the backward pass with a differentiable surrogate, such as a seoothed function or the identity.
Physics-informed approaches use automatic differentiation to compute derivatives of network outputs with respect imput coordinates. While automatic differentiation can handle non-differentiable components, they be problematic for innrrme- losses. Non-smoothness are undefined or discontinuous, and these issucs can propagate through the chain rule when higher-order derivalives are needed for PDE residuals. The result can be abrupt residual changes, numerical instability, and convergence to suboptimal solutions. Differentiable components are therefore preferable since they provide continuous derivatives, ensuring smooth residuals that better capture the physics. This reduces numerical artifacts, improves gradient quality, and supports stable optimization, leading faster convergence and more accurate solutions, with greater reliability.
create regions where derivatives
34


Reviews:
[{"review_id": "YmEUzGPYFA", "reviewer": "Reviewer_r2Xe", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. Mollifying the neighborhood kernel is a clean fix that makes GNO/GINO compatible with AD-based physics losses on arbitrary domains.\n2. There are clear improvements over finite differences for training with physics losses, strong point-cloud results, and a credible inverse-design demo with solver verification.\n3.  The ablations on radius/weights, subsampling of physics points, and notes on aggregation offer actionable tuning advice.\n\nWeaknesses:\n\n1. AD is exact for the computation graph, but gradients can still be affected by discretization, neighborhood changes at the support boundary, and non-smooth components. The phrasing should be carefully revised to be clear.\n\n2. Missing discussion on non-smooth operations: many operator stacks use ReLU/max/clipping and discrete neighbor updates. Claims about higher-order derivatives and smoothness need to address these cases explicitly. It might be that this is already taken care of by the Authors, by I could not understand this from the paper.\n\n3. Missing aggregation and radius sensitivity discussion: the paper touches on sum vs mean and shows a weight ablation, but a broader view of aggregation schemes and a clearer accuracy–cost frontier versus radius/neighbor count would help to improve its quality.\n\n4. The “Weight functions” paragraph drops in abruptly without re-stating purpose and design requirements.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes mollified Graph Neural Operators (mGNO) and a mollified GINO (mGINO) that replace the hard neighborhood indicator with smooth, compactly supported weights, enabling automatic differentiation (AD) of physics losses on irregular geometries and point clouds. This allows training with PDE residuals at arbitrary query points and supports inverse design by backpropagating through geometry. \n\nExperiments span 1D Burgers (on grids), nonlinear Poisson and hyperelasticity (on unstructured point clouds), and an airfoil optimization. Across tasks, AD with mGNO is reported to reduce training error relative to finite differences, outperform a meta-learning baseline on irregular domains at comparable runtime, and achieve notable speedups versus a numerical solver at similar accuracy; the airfoil design is validated with a high-fidelity solver.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I think that overall the paper is very interesting and well sectioned and written. There are some changes that need to be done before acceptance:\n\n1.  Clarify derivative claims, for example by replacing “exact derivatives/gradients” with phrasing such as “AD-exact with respect to the differentiable computation graph,” and briefly discuss discretization effects, support-boundary/neighbor changes, and their practical impact on stability and smoothness.\n\n2. Address non-smooth components. Clarify how you handle ReLU/max/clip and discrete neighbor updates: which derivative notion is used at kinks, whether smooth activations were considered for higher-order physics losses, and any observed differences (a small ablation would suffice).\n\n3. Differentiate from “physics-informed approaches without physics losses.” In that paragraph, explicitly state how your method differs from architectures that encode physics via constraints/projections/operator design rather than explicit PDE residuals. Note when each route is preferable (e.g., memory/compute vs flexibility on irregular domains).\n\n4. Weight-functions context. Before listing specific weights, add 1–2 sentences summarizing purpose (mollify hard neighborhoods for stable AD), key properties (smoothness, compact support, controllable locality), and default recommendations; note how radius is chosen (e.g., target neighbor count).\n\n5. For completeness, please add the initial conditions for Burgers equation solution in Figure 4.\n\n6. Regarding subsampling and data efficiency, please discuss recent training strategies such as:\n\n- Multiscale Training of Convolutional Neural Networks\n\n- Towards efficient training of graph neural networks: A multiscale approach\n\nThat also consider the data efficiency and time efficiency of training neural networks, with some overlap to the considerations in this paper.\n\n7. Since your method shows better approximation of derivatives, it would be interesting to quantify its stability or robustness. Can the authors perform an experiment with noisy inputs to show the performance of their method vs finite differences ?", "broader_impact_concerns": "I am not aware of negative impacts of this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "fCRQkoRowA", "reviewer": "Reviewer_eweb", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\n1. Clear, well‑motivated problem statement. The authors pinpoint a concrete bottleneck in existing neural operators: the binary neighbor indicator breaks automatic differentiation, preventing the use of gradient‑based physics losses and downstream design tasks.  By targeting this pain‑point the paper delivers a focused, easily understood contribution that resonates with both ML and computational‑physics audiences.\n\n2. Comprehensive experimental coverageThe evaluation spans:  (i) 1‑D viscous Burgers, (ii) 2‑D nonlinear Poisson, (iii) 3‑D hyper‑elasticity, and (iv) a real‑world NACA air‑foil inverse‑design task.  The breadth demonstrates robustness across equation families, dimensionalities, and mesh irregularities, while the air‑foil case illustrates real engineering impact.\n\n# Weakness\n\nIncremental methodological novelty. The core idea—smooth a discontinuous indicator—can be expressed in a single line of code.  Reviewers looking for deeper algorithmic insights or theory may find the contribution thin relative to the empirical section.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Enabling Automatic Differentiation with Mollified Graph Neural Operators introduces mGNO, a smoothed version of the Graph Neural Operator whose differentiable kernel weights let Autograd deliver exact spatial derivatives on irregular meshes and point clouds. Placing mGNO inside the Geometry‑Informed Neural Operator gives mGINO, enabling training with high‑resolution physics losses and little (or no) data while preserving differentiability across complex geometries. In experiments, Autograd‑driven mGINO cut Burgers‑equation data error by 20 × versus finite‑difference PINO, and on nonlinear Poisson and hyperelasticity tests, it achieved 2–3 orders‑of‑magnitude lower error while running 20~4 000 × faster than Meta‑PDE baselines and finite‑element solvers. The fully differentiable pipeline also supports inverse problems: optimizing an airfoil through mGINO produced a drag‑to‑lift ratio of 0.09, surpassing earlier neural‑operator designs. Overall, the paper shows that mollifying the GNO kernel unlocks exact gradients, substantial accuracy gains, major speed‑ups, and strong data efficiency for physics‑informed neural operators, extending their reach to complex geometries and design optimization tasks.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The single paragraph on future work (p. 17) is brief and not critical.\nIdeas like adaptive or hybrid losses appear tangential to the core contribution of smoothing the kernel weights in Eq. 6.\nA richer outlook, such as learning data‑driven mollifiers, coupling with adaptive meshing, or extending to multi‑scale operators—would better demonstrate the method’s potential and improve its novelty.", "broader_impact_concerns": "Beyond academic PDE benchmarks, the airfoil inverse‑design experiment underscores the method’s practical relevance.\nLeveraging fully differentiable aerodynamics, mGINO optimizes the wing profile, highlighting its potential for industrial‑scale design workflows that would otherwise depend on expensive CFD loops, matching or beating domain‑specific surrogate optimizers without additional heuristics.\nThis showcases a valuable new capability unlocked by the technique.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cfV2MYkYrJ", "reviewer": "Reviewer_VJVU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The paper addresses the increasingly discussed topic of combining data-driven and physics-driven approaches; this alignment with current research trends is noted.\n- Replacing the non-smooth indicator with a smooth weighting (mollifier) function is a standard idea from functional analysis and does offer a clean way to make the kernel differentiable.\n- The airfoil inverse-design task is well chosen and may interest both the CFD and machine-learning communities.\n\nWeaknesses:\n\n- Beyond substituting a predefined smooth weight for the indicator function, the implementation changes are minimal. The core algorithmic steps remain essentially the same as existing GNO formulations, limiting the incremental contribution.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors substitute the original non-differentiable indicator function in the GNO kernel integral with a smooth, differentiable weighting function. This refinement makes the entire GNO operator differentiable end-to-end, enabling automatic differentiation to compute exact spatial derivatives at arbitrary query points. As a result, physics-informed loss terms based on the governing PDE can be incorporated seamlessly during training.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- **Baseline selection is insufficient.** The comparison set omits several strong physics‐informed solvers. In particular, conventional PDE solvers such as Transolver [1] augmented with finite-difference PDE losses should be included as baselines; without them, it is hard to gauge the true benefit of the proposed method.\n- **Alternative avenues for incorporating equation information are not explored or discussed.** Beyond loss-based enforcement, recent approaches inject physics via conditional architectures (e.g., Unisolver [2], which conditions on all available equation components) or by casting the PDE itself as a computation graph (e.g., PDEformer [3]). A thorough empirical and conceptual comparison against such strategies is necessary for a fair assessment.\n- **Scalability to higher-dimensional, time-dependent problems remains unclear.** It is not evident how the current formulation extends to 2D + time PDEs—for instance, the classic Navier–Stokes dataset used in the original FNO paper [4], where a model ingests the first ten frames of a fluid sequence and predicts the next ten. Concrete guidance on applying mGNO in this setting is missing.\n- **Real-world validation is lacking.** Evaluating the model on at least one realistic, noisy dataset would help demonstrate robustness and practical value beyond synthetic benchmarks.\n\n[1] Transolver: A Fast Transformer Solver for PDEs on General Geometries\n\n[2] Unisolver: PDE-Conditional Transformers Are Universal PDE Solvers\n\n[3] PDEformer: Towards a Foundation Model for One-Dimensional Partial Differential Equations\n\n[4] Fourier Neural Operator for Parametric Partial Differential Equations\n\n###", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)