Neuro-Spectral Architectures for Causal Physics-Informed Networks
Authors: Arthur Bizzi
Leonardo M. Moreira
Márcio Marques
Leonardo Mendonça
Christian Júnior de Oliveira
Vitor Balestro
Lucas dos Santos Fernandez
Daniel Yukimura
Pavel Petrov
João M. Pereira
Tiago Novello
Lucas Nissenbaum
Venue: NeurIPS 2025 poster ()
Keywords: Physics-informed neural networks
spectral methods
causality
Primary area: machine_learning_for_sciences

Abstract: Physics-Informed Neural Networks (PINNs) have emerged as a powerful frame-
work for solving partial differential equations (PDEs). However, standard MLP-
based PINNs often fail to converge when dealing with complex initial value
problems, leading to solutions that violate causality and suffer from a spectral
bias towards low-frequency components. To address these issues, we introduce
NeuSA (Neuro-Spectral Architectures), a novel class of PINNs inspired by classi-
cal spectral methods, designed to solve linear and nonlinear PDEs with variable
coefficients. NeuSA learns a projection of the underlying PDE onto a spectral
basis, leading to a finite-dimensional representation of the dynamics which is then
integrated with an adapted Neural ODE (NODE). This allows us to overcome
spectral bias, by leveraging the high-frequency components enabled by the spectral
representation; to enforce causality, by inheriting the causal structure of NODEs,
and to start training near the target solution, by means of an initialization scheme
based on classical methods. We validate NeuSA on canonical benchmarks for lin-
ear and nonlinear wave equations, demonstrating strong performance as compared
to other architectures, with faster convergence, improved temporal consistency
and superior predictive accuracy. Code and pretrained models are available in
https://github.com/arthur-bizzi/neusa.

Neuro-Spectral Architectures for Causal Physics-Informed Networks
Arthur Bizzi', Leonardo Moreira 3 Márcio Marques', Leonardo Mendonça', Christian Oliveira", Vitor Balestro", Lucas Fernandez', Daniel Yukimura², Pavel Petrov?, João M. Pereira Tiago Novello?. Lucas Nissenbaum
ficole Polytechnique Fédérale de Lausanne (EPFI.). Instituto de Matemática Pura Aplicada (IMPA), 3 Universidade do Estado do Rio de Janeiro (UERJ), "Laboratório Nacional de Computação Cientifica (LNCC), University Georgia (UGA)
Abstract
Physics-Informer Neural Networks (PINNS) have emerged a powerful frame- work for solving partial differential equations (PDEs). However, standard MLP- based PINNs often fail to converge when dealing with complex initial value problems, leading to solutions that violate causality and suffer from a spectral bias towards low-frequency components. To address these issues, we introduce NeuSA (Neuro-Spectra) Architectures), novel class of PINNs inspired by classi- cal spectral methods, designed solve linear and nonlinear PDEs with variable coefficients. NeuSA learns projection of the underlying PDE onto spectral basis, leading to a finite-dimensional representation of the dynamics which then integrated with adapted Neural ODE (NODE). This allows us overcome spectral bias, by leveraging the high-frequency components enabled by the spectral representation; enforce causality, by inheriting the causal structure of NODEs, and to start training near the target solution, by means of an initialization scheme based on classical methods. We validate NeuSA canonical benchmarks for lin- ear and nonlinear wave equations, demonstrating strong performance as compared other architectures, with faster convergence, improved temporal consistency and superior predictive accuracy. Code and pretrained models are available httpss://github. com/ arthur- bizzi/neusa
Introduction
The introduction of Physics-Informed Neural Networks (PINNs) has sparked interest in using neural networks solve partial differential equations (PDEs) 2 PINNs enable data-efficient modeling of complex systems by embedding physical laws directly into the loss landscape. This approach has opened new possibilities in scientific computing, with applications spanning wide range of subjects, including fluid dynamics and climate modeling [5 61. biomedical simulations
3. material science (9) 10 M. and others [12]
13
14
15
16
We consider the following
initial-boundary value problem for E [] and x E Ω
C RR::
n(t, x) = F (t, VU, VVn). u(0,x) = ux(x)
dt
where [0, x Ω R denotes the (vector-valued) solution, Vu and V denote its first- and second-order spatial derivatives, and smooth function. un(x) is the initial condition (Cauchy data), and assume some boundary conditions are imposed OΩ. PINNs consist of
"Corresponding author: artur.continhohizzi@eep...h
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Refative Error
GT
PMS
Ores
FLS
PiNksFomer
Ne.5A lours)
P/M Onic
- 0.0
107
Aeuus Souind
-1.00
1.0
0.000
-0.24
101
2.0
0.00
IMSE
0.5469
0.1978
0,3885 1.4231
0.1296
0.32
10:1
250
500
750
1000
1250 1500
Training time
Figure I: We present NeuSA, a theoretically grounded Physics-informed neural architecture. On the left, we compare various models on wave propagation problem. The dashed lines represent the discontinuities a stratified heterogeneous medium. NeuSA achieves the lowest relative error (rMSE) and most accurately preserves sharp wavefronts and reflections. On the right, we show the evolution of the relative L2 error during training. NeuSA converges more rapidly and consistently.
0 ppaaametrized coordinate networks ug [0, T] 2 R", trained to approximate the solution by minimizing a composite physics-informed informed loss function ((0). This loss typically includes terms for the PDE well as for the initial/boundary data, e.g..
d (PDE(0)= = di us (t ug, Vug, VVna) Lic(0)= - Arc /luolllu0,00000 x) - (x)] (2)
where || ||2 denotes the regular norm in 1.2 ( 7] x Q), and Aic is the weight for the associated residue term (a term LBc representing boundary conditions could be also included). PINNs thus provide a flexible, data-driven, and meshless framework for approximating PDE solutions using neural networks.
Data-Driven. As neural networks, PINNs are particularly suited for handling heterogeneous, noisy, or incomplete measurement data, which can be efficiently combined with physical priors via physics- informed losses [4 13 17
18 19
As coordinate networks, PINNs represent continuous interpolations of the underlying solutions and may be evaluated at arbitrary spatial or temporal coordinates. In higher- dimensional problems, this property can be combined with random sampling strategies to substantially reduce the number of samples needed to approximate the solution 20
Mesh-Independee
However, standard PINNs often struggle to enforce fundamental structural aspects of the underlying solutions. In most cases, they rely on general-purpost feed-forward architectures, such as the standard Multi-Layer Perceptron (MLP) [21]. or on specialized MLP- based variants that enhance expressivity through indivaatio-functtioo modifications, including QRes (22). FLS (23]. This generic structure design often leads to issues related to spectral bias, causality and limited generulization capacity:
Spectral Bias. Regular coordinate networks based sigmoid or rectifier activations often struggle represent high-frequency components, leading to issues with representing detailed and/or multi-scale solutions (24 25). This effect is often mitigated with Fourier-Featur (FF) layers, sinusoidal encoder layers designed to inject high-frequency representations into a network's 's architecture [26). Still, FF layers require fine-tuning to avoid overfitting and noise.
Causality. PINNs are notorious for violating causality and temporal consistency due to their simultaneous training over the entire time domain [17]. These issues may manifest in the form of incorrect initial conditions or non- non-physical convergence to trivial solutions. Attempts have also been
made minimize these effects with modified losses [17 [9]
Generalization Capacity. bl.ed PINNs may struggle with extrapolation beyond their training domain 27 28), which has been tackled with alternative training strategies (29].

---

Due to these issues, PINNs often fail to converge to the true solution when solving complex time- dependent problems. Instead, they may overfit and converge trivial equilibrium solutions. Such shortcomings are common when solving problems with strong time dependence, as evidenced by their relative lack of success when applied to linear and nonlinear wave equations [12 31 This
stands in stark contrast to PINNs capacity for solving parabolic and elliptic equations
32
We propose NeuSA a new family of Neuro-Spectral Architectures designed for solving space- inhomogencous and/or nonlinear time-dependent PDEs. NeuSA uses the spectral decomposition to obtain a method-of- lines 33 34 discretization of PDE into large system of ODEs, which is then modeled using Neural ODE (NODE) 335 (see Figure 2). Figure showcases NeuSA 's results on 2D wave propagation task, demonstrating significant improvements over prior methods in accuracy, speed, and temporal consistency. Our contributions may be summarized as follows.
Causality. NeuSA is apectral-method -assed architecture for neural PDEs, such as [36]. Consequently, inherits the causal structure of classical methods, including exact initial conditions and uniqueness, while retaining data-friendly, mesh-less representation.
Spectral fidelity. The choice of global spectral bases allows NeuSA to overcome the spectral bias commonly attributed to MLP-PINNs, offering a Fourier-Feature Networks.
alternative to
Analytical initialization. The interpretable structure NeuSA as a neural extension of spectral methods enables specialized initialization schemes in which networks are initialized as the solution of closely related linear homogeneous problems, at no training cost.
Time-extrapole extepoaxrapoottion. Due to its causal formulation, NeuSA displays strong time-extrapolating performance, enabling simulation beyond training intervals.
((0)
a(h))
((tz)
u(12, x)
bix)
Figure 2: Neuro-Spectrr Architecture. Above: The spectral coefficients for the initial conditions i(0), flowing according to NODE. Below: The spatial input being encoded into the spectral basis functions b(x). Coefficients and bases are then combined to yield the final result.
1.1 Related work
Several works have explored designs for physics-informed neural machine learning, including large body of literature on operator learning 337 28 Multiple recent works on Physics- Networks propose alternative architectures for representing solutions, enhancing their expressiveness, spectral representation, or temporal coherence.
Quadratic Residual networks (QRes) [2) introduce class of parameter-efficient neural networks by incorporating a quadratic term into the weighted sum of inputs before applying the activation functions each layer of the network. This modification enables QRes to approximate polynomials using shallower networks, resulting in compact yet expressive models.
Layst-Layer Sine (FLS) [23, also referred to as sf-PIN, introduces sinusoidal encoding layers within the PINN framework. FL. nets utilize sinusoidal encoding layers, in an approach closely analogous to that of Fourier-feature networks, to mitigate spectral bias and enhance input gradient distribution.
3

---

PINNsFormer adapts the transformer architecture to the PINN setting, leveraging attention mechanisms to model temporal dependencies among state tokens, with the goal of achieving enhanced temporal consistency.
Neural Ordinary Differential Equations (NODEs) 35 class of 'continuous-depth" residual networks that model inference as the integration of a continuous-time process, effectively solving an ODE whose dynamics are parameterized by neural network. NODEs have proven powerful for modeling continuous-time dynamics and have been applied to physics-informed learning, gen-
erative modeling, time-series forecasting, and morphing 38 39
43 However, their highly sequential numerical structure makes them equivalent to ultra-deep residual networks, which can significantly slow down training. As a result, their application PDEs remains relatively underexplored [44]
2 Neuro-spectral architectures
Neuro-Spectral architectures are defined as models that employ spectral decomposition to reduce a PDE defined over dimeiieddimeniioall space to ODE system in finite dimensions, subsequently training a NODE to approximate the latter using a physics-in informed loss. This formulation interprets I abstract Cauchy problem over the Hilbert space L2())), treating the solution u(t,-) as a time-paramen pareeearraterzee family u : IR >> (2)))) The spectral decomposition [45] 46 47 consists of approximating u(x. f) by its projection onto the subspace spanned by finite subset b(x) of an orthonormal basis of L (22)". Given a truncated spectral representation with harmonics C1. C2,.... Cd.
the solution u can be expressed in terms of expansion over elements of the basis tensor b(x):
Ω - Cc. Xce, whose coefficients form tensor ú(l) [0, T] + C X-xxxx leading to
Σ (x)b))x( (x(()) = / u(t, ())x(x)xx
u(t,x) x)
Where lt denotes a dddimeensional multi-index, (4) : [0, 7 - Rm and by (x): : + R represent the k th indexed element in ú and b, respectively. Substituting (3) into (1 leads to method-of-lines discretization, resulting in an ordinary differential equation for the coefficients:
33
ú F(ú).
dt
where F: CC (CexX C xxxxx usually does not admit simple closed-form expression
for a general basis b. Instead of deriving it explicitly, we learn F as parameterized network Fg. Inference in Neuro-Spectral model then proceeds as follows (see(Figure 3 and [48] ):
1. Project the initial conditions onto an orthonormal basis. Sample the initial conditions u(0, x) densely and extract their spectral representation ú in terms of the basis b: u(0, x) = Dn úg (0)bb)))))
2. Integrate coefficients in time according to NODE. Use the coefficients tensor ú as input to
NODE with vector field Fg and integrate with a high-order method: ún (())+ Fo (ū(T))
3. Reconstruct the solution and perform training. Multiply the obtained coefficients úa by their corresponding basis functions (x) to obtain u(t, x). This representation can be differentiated analytically to compute physics- inforredd losses: un (f, x) = Σx (x))
2.1 Spectral decomposition and initialization
The choice of the basis b can enforce specific properties of the solution, as well as ensure the fulfillment of the given boundary conditions. In this work, we initialize b as the Fourier basis and its odd and even extensions in terms of the sine and cosine functions. This choice enables the representation of homogeneous periodic, Dirichlet, and Neumann boundary conditions for rectangular domains. The spectral projection ú can then be computed in several ways, depending on how the spatial domain is sampled (see Appendix B). We adopt the Fourier basis for two main reasons: first, to overcome spectral bias, inspired by the success of Fourier-Feaurr layers; and second, to allow for the simple and accurate representation of linear translation-invarian (LTT) differential operators as scalar multipliers [49]. We use the latter to implement an improved initialization scheme for the NODE, described as follows.

---

((())
GG4)
Figure 3: Inference in Neuro-Spectral model. The initial conditions are decomposed into their spectral coefficients, which are propagated time via a NODE. The time-iterated coefficients are then reconstructed into the solution at later times.
1. Linearize the PDE. Extract linear translation invarriant approximation for F:
d u = (u, Vu, VVu.... aou + Σ
d²
u
+
Σ
(5)
dt
dx;
dx;dx;
1.
Fourier multiplier. Derive the associated Fourier multiplier M C Xe2xxxxxx defined as an element-wise polynomial on the k-th Fourier frequency corresponding to the k-th harmonic:
d
d dt
F Lincar
M
(6)
where O stands for the Hadamard (element-wise product.
3. Initialize the vector field F,. Initialize NODE near this approximate linear solution by augmenting the learned vector field with the analytical multiplier:
Fú = (M @ û) + (F))
where (ū) is a neural network initialized with mean zero and unit variance, and small parameter. In this way. the network starts close to the solution of the associated LTI problem, serving a strong prior for training. During optimization, Fe learns compact representation for the dynamics, effectively leading to a neural generalization of
non-linear and/or aun-translation-invee the classical spectral method. See Section for explicit examples, and Appendix B.2 for the detailed architecture of FB.
2.2 Neural ODE, time integration and causality
As discussed, neuro-spectral models rely NODE propagate the spectral coefficients forward in time. The vector field to be integrated over as part of inference composed of the cear-analytical initialization discussed above along with multilayer perceptron Fs:
d ú = (4)
o + eF)
dt
The NODE receives as input a tensor of size C1 X C2 X ... Ca n. corresponding the first time-slice of the solution, and outputs Exmmlees slices in the form of tensor of shape samples x 22 ca X n. This special treatment of the time dimension, characterized by the sequential nature of integration, is what equips NeuSA with causal structure. In fact, it is possible to prove NeuSAs are flows [50 51 as summarized in the following theorem:
Theorem 1. For band- limited initial conditions un and globally Lipschitz neural vector fields Fo. the orbits created by NeuSA satisfy the initial conditions and uniqueness:
1. fulfillment of initial conditions: (0, x) u(0,x);
2. uniqueness: (0, -) + 3 (0.-) = (t, (t,-t)) Vt E [0,
A more detailed exposition of this theorem as well as its proof may be found in Appendix A It is. in essence, result of the properties of flow operators for ODEs combined with the uniqueness of the spectral decomposition for band-limiter functions. This holds by construction, regardless of training.
5

---

2.3 Losses and training
Training NeuSA similar to training a common MLP-PINN. The main difference is that we can no longer differentiate directly with respect to time, as it is no longer an input coordinate; it is instead implicitly encoded as the time-steps for the NODE iteration. Nevertheless, time and space derivatives may be calculated in straightforward manner:
d dx,
d
di uo x)
((o))(((()))(((
(o,
Σ Úg (1) dx,
bx(x),
where by construction ((( (t) d (t). Note that the cost of calculating derivatives does not increase meaningfully for higher-order spatial derivatives, as opposed to the exponential increase in computational cost incurred by naively stacking derivatives with autograd [52]. We may then sample the domain 52 and evaluate the associated Physics-Infiee Informed residue with
at
CPDE (0) = -E Σ 6EJO,2 XyEn
us ua xj) - F(4,, Xy, ug Vus, VVus) di
d
(10)
where i denotes the i-th integration time step for the NODE and xj denotes the coordinates at the j-th spatial sample point. Note that NeuSA automatically complies with initial and boundary conditions and therefore does not require loss terms for them.
Note that neither time nor space samples need to be uniformly distributed; nevertheless, space samples must remain constant across all times for each pass, as opposed to conventional PINNs. This comes the advantage that a single forward pass necessary to evaluate the loss over all samples, as opposed to the multiple passes needed for common PINNs. This will allow NeuSA architectures to achieve training speeds comparable to those of purely neural approaches, despite their reliance on computationally intensive NODE integration.
3 Experiments
We evaluate NeuSA on boundary and initial value problems for three PDEs: the 2D wave equation, the 2D Burgers' equation, and the 1D nonlinear sine-Gordon equation. In all cases, we address the forward (direct) problem, where the models are trained to learn approximate solutions given known conditions. We compare the performance and accuracy of NeuSA against several established MLP-based PINN architectures: the original PINN M. QRes
FLS (33, and PINNsFormer
Training setup. NeuSA and all baseline models are implemented in PyTorch 53 The baseline configurations follow the setups described in PINNsFormer 4 and RoPINN K PINN, QRes, and FLS are initialized using Xavier initialization [54), with the hyperbolic tangent as the activation function (except for the first FLS layer, which acts as a Fourier feature mapping [55] 26)). All remaining hyperparameters were tuned to achieve the best performance for each model (c.g. weights for the initial and boundary condition losses). For the NODEs used NeuSA, we adopt the implementation given by the TorchDyn library (56). The vector fields Fg modeled as MLPs with dimensionwise layers (see Appendix [B.2) with two hidden layers, ReLU activations, and Glorot initialization, and are integrated using fourth-order Runge Kutta solver. Gradients computed via standard backpropagation through the ODE solver.
Training is performed using the Adam optimizer [57]. NeuSA's strong architectural priors enable the use of larger learning rates compared to the baseline models, which are trained with the recommended rate of 10 (see 1 of (58). All bascline models exhibited reduced performance when trained with learning rate 10 2 due to instability. Each experiment was run several times, and the mean of each evaluation metric
reported. To evaluate model accuracy, we consider the standard rMSE (relative mean squared error, also called the relative L2 error) and rMAE (relative mean abxolute error, also called the relative L1 error) metrics. The ground-truth solution for each problem is obtained with high accuracy numerical solvers (see (Appenddi C). Experiments were executed on an Nvidia RTX 4090 GPU (24 GB VRAM). Results are summarized in Table
Additional details and extended results for each experiment are provided in Appendii C and additional experiments exploring NeuSA's training and inference time may be found Appendix D.

---

3.1 Wave equation: 2D Layers, 3D Layers and Marmousi
As our first benchmark, we consider an initial-value problem for the linear wave equation in an infinite heterogeneous medium, canonical problem in acoustics, seismics, and electromagnetism, [59 60) 61]:
0 d12 u = (x) u(x,0) exp
|x|2 2x
a u(x, 0)
(11)
where A denotes the Laplacian in the spatial dimensions, = 0.1, and denotes the material- dependent wave speed. This material heterogeneity leads to wave reflection and refraction at the interfaces between layers. We evaluate NeuSA in three scenarios of increasing complexity: the 2D wave equation for an environment with three horizontal layers (hereafter referred to as 2D Layers), the 3D wave equation in medium with two horizontal layers (here referred to as 3D Layers), and the 2D Marmousi model. all cases, the spatial domain is truncated to [-2, 2] for = 2, 3. This is valid approximation to the infinite domain scenario long the propagating waves do not reach the domain boundaries.
For the 2D Layers and 3D Layers cases, the propagation medium consists of three and two horizontal layers, respectively. Such simplified test problems are frequently found in seismic datasets [62]. The Marmousi reservoir model 663 is canonical benchmark with highly complex stratified medium in two dimensions, containing folds rocks and an overlying water layer (64) Images depicting both 2D media may be found in Appendix C
For NeuSA, we initialize the neural vector field near the solution of the homogeneous wave equation and constrain Fg to take a low-rank form, as detailed in Appendix B The model is trained for 2,000 steps, using 201 basis elements and integrated over 201 time steps. The baseline models are trained
for 20,000 steps, using 10,000 random collocation points for the PDE and 1,000 points for the initial
conditions. See FFiuuee for an image depicting the results obtained for the 2D Layers and Figure for the Marmousi reservoir. The results obtained for the 3D Layers case are presented in Appendix C
Ground-Truth
PINN
FLS
QRes
NeuSA
0.66
6.0
Amplitude
t= 11133
0.00
-0000
t= 2.00
rMSE
0.6310
0.7186
0.4451
0.1355
Figure 4: Results for the wave equation over the Marmousi benchmark. NeuSA is able to achieve solution that is much closer than reference methods the ground truth.

---

3.2 Sine-Gordon equation
In this example, we consider the 1D sine-Gordon (s-G) equation, a generalized nonlinear wave equation that has applications in soliton collisions and inverse scattering Consider the problem
Ou du 10))nnu), θt
1 exp V2TT
(-
|x| a 20
x(0,x)
u(x,
(12)
dz
ōt
with (x, E [-4,4] [0, 3). zero-Dirichlet boundary conditions at =4, and 0.1.
For the numerical experiments, NeuSA is trained with 201 frequencies and 201 time steps, using 1, 000 steps with learning rate of 0.01. The PINN, QRes, and FLS models are trained on regular grid of dimension 201 x 201, while 101 x 101 grid considered for the PINNsFormer model. All models other than NeuSA were trained for 10, 000 Adam steps. The computational results are visualized Figure
Ground buth
10
Paztstomer
R.S
0.5
0.-
PINN
0.0
Offirs
NeuSA
-0.5
0.5
-1.0-
0.05
0.35
0.30
0.00
0.05
-001
0.05
0.00
-002-
-0.10
-0.05
-003
-004
-0.30
Figure 5: Results from the sine-Gordon equation. Above: the solution versus the ground truth. Below: the residues between them. NeuSA both faster and more accurate than the baselines.
3.3 2D Burgers' equation
We address 2D generalization of Burgers' equation, commonly used benchmark in computa- tional fluid dynamics (66) to model the development of two-dimensional shocks [67]. We consider the initial-boundary value problem
a u - Vu + vdu, 01
u(x, 0) (sin(TT)) (T))) sin(TTr2, cos(TTI)) cos(T..2))
(13)
where (u(x, t), u(x, t)), (x, t) [0, 4 x [0, 1). and 0.01, and periodic boundary conditions in x are assumed.
NeuSA is trained for 200 steps using 201 x 201 components over 201 time steps, with vector field initialized near the solution of the corresponding heat equation. The baseline models are trained for 20,000 steps, using 10,000 collocation points for the PDE, 1,000 for the initial condition, and 500 for the boundary condition.
We conducted an additional experiment on Burgers' equation to assess our method's ability to extrapolate the approximated solution time beyond its training interval. All models were trained on the interval [0. 1] and then evaluated on the extended interval [0. 2]. We analyzed performance in successive time instants, quantifying how prediction accuracy degrades as the time gets further from the training region. Because our method integrates learned vector field over time, we expected it to extrapolate smoothly into future time, rather than exhibit the localized overfitting often observed in standard MLP- based architectures. The results shown in Figure 6 confirm expectations.

---

Ground Truth 1-2.0
QRes t-2.00
Relative L2 Time
0.50
3
0.25
144
NeuSA
0.00
1.2
QRes
Rolative Frtor
-0.50
Training region
Extrapolation region
0.8
NeuSA t=2.00
0.2
PINN V: t=2.00
0.4
0.6
0.1
0.2
22
22
0.0
0.2
-0.1
0.0
3
0.00
0.25
0.50
0.75
1.00
3.25
1.50
1.75
2.00
Figure 6: Extrapolation results for the Burgers' equation. Left: solutions from PINN and QRes degenerate strongly when extrapolating; in contrast, NeuSA's solution retains qualitative behavior observed the Ground Truth. Right: NeuSA outperforms all the baselines when extrapolating, maintaining comparable performance. Other well performing models, such QRes, quickly diverge.
Table I: Average experiment metrics for each model. TT refers to the training time in seconds. Despite
being generally trained for less time, NeuSA consistently matches outperforms the baselines.
2D Layers
Marmousi
3D Layers
2D Burgers
Sine Gordon
Model
MMSE
TT
rMSE
TT
rMSE
TT
rMSE
TT
rMSE
TT
PINN
0.545
566
0.698
635
0.073
5990
0.221
871
0.139
976
QRes
0.115
750
0.412
718
0.021
8775
0.073
135
0.020
1,315
FLS
0.590
577
0.684
648
0.070
6179
0.202
885
0.135
1,015
PINN&Former
1.072
1484
1.053
2, 294
0.681
3,
NeuSA
0.075
530
0.171
573
0.008
702
0.051
112
0.001
215
Results and analysis
The results for the wave equation show that NeuSA can accurately propagate waves in complex media, leading to an error between one and two orders of magnitude smaller than the baselines. This is despite NeuSA being trained for 10 times fewer training steps. It is worth noting that key challenge in seismic and acoustic simulations with heterogeneous media is to capture complex wave patterns generated by reflections at material discontinuities. NeuSA is the only approach that accurately recovers second-order reflections, as shown Figure at
The results for the sine-Gordon equation in Figure 5) iidiicate that NeuSA very quickly learns the dynamics for the system with outstanding precision, obtaining error metries an order of magnitude smaller than QRes, the next best result, with a substantially smaller training time. NeuSA's causal structure allows the ground-truit truth. In contrast, PINN, FLS and PINNsformer fail to do so, despite being trained for 10 times more epochs.
propagate the solution from initial conditions, staying in close agreement with
The results for the Burgers' equation show that NeuSA successfully captures the qualitative behavior of the solution, including the formation of detailed two-dimensional shocks. Quantitative results in Table show that NeuSA exhibits strong performance, while requiring substantially less training time, leading gains over PINNs, FLS, PINNsFormer and QRes.
These results are underscored by the very strong performance in the extrapolation task shown in Figure 6] which demonstrates that NeuSA has learned robust internal representation for the system's autonomous dynamics. As a result, NeuSA can estimate and extrapolate solutions over large times intervals beyond the training domain, leading to substantially superior performance compared to the baselines. We attribute this capability to NeuSA's causal structure.
Although built on the inherently slower NODE framework, NeuSA trains significantly faster than MI.P- PINNs. This efficiency likely due to the physical and causal priors built into the architecture, allowing for convergence in considerably fewer steps. Moreover, neuro-spectral models

---

compute the solution across the entire domain in single forward pass, in contrast to the thousands of passes needed for PINN, reducing the relative computational overhead. [ppendix D has longer discussion on this effect. Additionally, NeuSA produces remarkably accurate solutions, which we attribute to its spectral representation capabilities and its automatic compliance with initial and boundary conditions. NeuSA does not need to train on boundary data, allowing to perform optimization only on the physics-informed loss. This helps avoid the conflicting gradients of data and equation losses, which often lead to unstable training, as shown in the rMSE-versus-time plots of the tested classical architectures in Figure
Nevertheless, the method does have limitations, which we discuss in greater detail in Appenddxx E First, all experiments take place on simple spatial domains, which allows using the framework of the Fourier bases; more complex geometries would require the implementation of more generic bases, which could hinder the initialization process. Second, NeuSA integrates the vector-valued solution ixxx f) using Munge-Kutta methods which, despite their efficiency, may become unstable when applied stiff problems. Finally, NeuSA's performance seems to be strongly influenced by its analytical initialization procedure; while this approach enables exceptionally fast training when suitable linear approximation is available, its effectiveness significantly diminishes when initialized with prior dynamical information.
5 Conclusion
We have introduced NeuSA (Neuro-Spectral Architectures), novel class of PINNS grounded in numerical spectral methods and designed to solve time-dependent nonlinear PDEs in inhomogeneous media. These architectures overcome the spectral bias and causality issues associated with MLP- PINNs by construction. Furthermore, we have shown that the architecture inherently satisfies boundary and initial conditions, as well as ensures solution uniqueness, without incurring any additional training cost.
Based the numerical experiments for the 2D wave equation in heterogeneous media, 2D Burgers' equation, and ID nonlinear sine-Gordon equation, we observed that NeuSA can solve complex problems, often achieving smaller error with significantly shorter training time than reference physics- informed architectures. Such capabilities are likely due to the improved initialization procedure and causal structure of NeuSA, which results in more physically relevant solutions.
The development NeuSA paves the way for future research and applications. In particular, the architectures be adapted to solve other challenging PDEs, such the widely studied Navier- Stokes equations [68). or even pseudodifffeential equations arising in wave propagation [69 70) Although initializing with Fourier basis has proven highly effective for the applications considered in this work, other bases could also be explored to enhance performance and, thus, make NeuSA less dependent on the mesh. Finally, other numerical methods could be used replace the fourth- order Runge-Kutta Kutta and improve the overall performance of NeuSA, especially when applied to stiff problems.
Acknowledgments
This work was supported by Petrobras. João Pereira is thankful for a start-up grant from the University of Georgia. We also acknowledge financial support from Google. We would also like to thank Deborah Oliveira and Lucas Schwengber for fruitful discussions.
References
[1] Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics-informee neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
[2] Gregory Kang Ruey Lau, Apivich Hemachandra, See-Kiong Ng, and Bryan Kian Hsiang Low. PINNACLE: PINN adaptive collocation and experimental points selection. In The Twelfth International Conference on Learning Representations, 2024.
[3] Haixu Wu, Huakun Luo, Yuezhou Ma, Jianmin Wang, and Mingsheng Long. RoPINN: Re- gion optimized panysics-iffommet neural networks. In A. Globerson, L. Mackey, D. Bel-
10

---

grave, A. Fan, U. Paquet, J. Tomezak, and C. Zhang. editors. Advances in Neural Information Processing Systems. volume 37, pages 110494-110532. Curran Associates.
Inc., 2024. wtw.hhtpss:/// p/tpceeed///oooeiingggg.s neur ips. wwwpappr.r... fiw.fflllf/plll pttee//pppper/2///1leee cr555aal000488000444445P2----ppp
pdff
[4] Zhiyuan Zhao, Xucying Ding, and B Aditya Prakash. PINNsFormer: A transformer-based framework for physics-informed neural networks. arXiv preprint arXiv:2307.11833, 2023.
[5] James Donnelly, Alireza Daneshkhah, and Soroush Abolfathi. Physics-informed neural networks surrogate models of hydrodynamic simulators. Science of The Total Environment, 912:168814, 2024. doi: https://doi.org// http////ooorr////// 10. 11.1000.1100
[6] Constanza A. Molina Catricheo, Fabrice Lambert, Julien Salomon, and Elwin van 't Wout. Mod- eling global surface dust deposition using physics-informed neural networks. Communications Earth & Environment, 5(1):778, 2024. doi: 10110000/43244-022--199494
[7] Stefano Buoso, Thomas Joyce, and Sebastian Kozerke. Personalising left-ventricular biophysical models of the heart using parametric physics-informed neural networks. Medical Image Analysis, 71:102066, 2021. doi: 10.1016/j www..e.......110002.
[8] Clara Herrero Martin, Alon Oved, Rasheda A. Chowdhury, Elisabeth Ullmann, Nicholas Pe- ters, Anil A. Bharath, and Marta Varela. EP-PINNs: Cardiac electrophysiology characterisation using physics-informer neural networks. Frontiers in Cardiovascular Medicine, Volume 8 2021, 2022. doi: 23338//few..2221.768888
[9] Shahed Rezaei, Ahmad Moeineddin, and Ali Harandi. Learning solutions of thermodynamics- based nonlinear constitutive material models using physics-informed neural networks. Compu- tational Mechanics, 74(2):333-366 366, 2024.
[10] Enrui Zhang, Ming Dao, George Em Karniadakis, and Subra Suresh. Analyses of internal structures and defects in materials using physics-informed neural networks. Science Advances, 8(7)eabk00444, 2022. doi: 10.1 26/sciadv.abk0644
[11] Bin Zheng, Tongchun Li, Huijun Qi, Lingang Gao, Xiaoqing Liu, and Li Yuan. Physics- informed machine learning model for computational fracture of quasi-britten materials without labelled data. International Journal of Mechanical Sciences, 223: 2022. doi: 10.1016/j.
ijmecsci.2022 1022.0
[12] Tim De Ryck, Siddhartha Mishra, and Roberto Molinaro. Weak physics informed neural networks for approximating entropy solutions of hyperbolic conservation laws. In Seminar für Angewandte Mathematik, Eidgenössische Technische Hochschule, Zürich, Switzerland, Rep, volume 35, page 2022, 2022.
[13] Ali Hasan, João M. Pereira, Robert Ravier, Sina Farsiu, and Vahid Tarokh. Learning partial differential equations from data using neural networks. In 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3962-3966, 2020. doi: 10.1109/ ICASSP40776.2022 217.662222.900.5500
[14] Zhiping Mao, Ameya D Jagtap, and George Em Karniadakis. Physics-informed neural networks for high-speed flows. Computer Methods in Applied Mechanics and Engineering, 360: 1 12789, 2020.
[15) Ravi G Patel, Indu Manickam, Mitchell A. Trask, Nathaniel A.and Wood, Myoungkyu Lee, Ignacio Tomas, and Eric C Cyr. Thernodynamically consistent physics-informed neural networks for hyperbolic systems. Journal of Computational Physicx, 449: 10754, 2022.
[16] Chengping Rao, Hao Sun, and Yang Liu. itysiie-informed deep learning for computational clastodynamics without labeled data. Journal of Engineering Mechanics, 147(8)004021033, 2021.
[17] Sifan Wang, Shyam Sankaran, and Paris Perdikaris. Respecting causality for training physics- informed neural networks. Computer Methods in Applied Mechanics and Engineering, 421:
1 116813, 2024.
[18] Liu Yang, Xuhui Meng, and George Em Karniadakis. B-PINNs: neural networks for forward and inverse pde problems with noisy data. Journal of Computational Physics, 425:109913, 2021.
Bayesian physics-informed
11

---

[19] Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. informed neural networks for forward and inverse pde problems. Computer Methods in Applied Mechanics and Engineering, 393:1 14823, 2022.
physics-
[20] Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non- adaptive and residual-based -aaee adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 116711, 2023.
[21] Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli. Scientific machine learning through -yysiis-inffrmed neural networks:
Where we are and what's s next. Journal of Scientific Computing, 92(3):88, 2022.
[22] Jie Bu and Anuj Karpatne. Quadratic Residual Networks: A New Class of Neural Networks for Solving Forward and Inverse Problems in Physics Involving PDEs, pages 675-683. doi:
10.1137/1.978161
[23] Jian Cheng Wong, Chin Chun Ooi, Abhishek Gupta, and Yew-Soon Ong. Learning in sinusoidal spaces with physics-informed neural networks. IEEE Transactions on Artificial Intelligence, (3):985-1000, 2024. doi: 10. 110//TA..022.39933566
[24] Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449: 10768, 2022.
[25] Zhi-Qin John Xu, Yaoyu Zhang, Tao Luo, Yanyang Xiao, and Zheng Ma. Frequency principle: Fourier analysis sheds light on deep neural networks. arXiv preprint 1rXvv2990..0653,, 2019,
[26] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barro,, and Ren Ng. Fourier features let networks learn high frequency functions low dimensional domains. Advances in neural information processing systems, 33:7537-7547, 2020.
(27] Keyulu Xu, Mozhi Zhang, Jingling Li, Simon S Du, Ken-ichi Kawarabayashi, and Stefanie Jegelka. How neural networks extrapolate: from feedforward to graph neural networks. In International Conference on Learning Representations (ICLR), 2021.
[28] Min Zhu, Handi Zhang, Anran Jiao, George Em Karniadakis, and Lu Lu. Reliable extrapolation of deep neural operators informed by physics or sparse observations. Computer Methods in Applied Mechanics and Engineering. 412: 11064, 2023. ISSN 0045-7825. doi: https:
//doi.org/10. 10. 1016 wmt.j.ccaaa...aaamm.a... 2023. 16624. URL https://s:///ww sciencedirect.. wwteeeetttttteeeeeeeeeeeeeeeeeee commsciennce////
article/pii/ Pttzicle/p3///04477555001113 //////000577222001883383
[29] Jungeun Kim, Kookjin Lee, Dongeun Lee, Sheo Yon Jhin, and Noseong Park. Dpm: A novel training method for physics-informed neural networks in extrapolation. Proceedings of the AAAI Conference on Artificial Intelligence, 35(9):8144) 8154, May 2021. doi: 10. 1609/aai... 13599.16992. Uttps:ss:://oos aaai org/ index. www.inppppppppppppeeeeppppppppte AAAI/ htt///rricleeeiiieee wtwp////ee/////////2//
[30] Yi Ding, Su Chen, Hiroe Miyake, and Xiaojun Li. mhysics-informed neural networks with Fourier features for seismic wavefield simulation in time-domale nonsmooth complex media. arXiv preprint arXiv:2409.0355 2024.
[31] Ben Moseley, Andrew Markham, and Tarje Nessen-Meyer. Solving the wave equation with thysics-informet deep learning. arXiv preprint arXiv:2006. 11894 2020.
[32] Makoto Takamoto, Timothy Praditia, Raphael Leiteritz. Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Information Processing Systems, 35:1596-1611 2022.
[33] Randall J. Le Veque. Finite Differrence Methods for Ordinary and Partial Differential Equations: Steady-State and Time- -eeendeendddtt Problems. Society for Industrial and Applied Mathematies. 2007. doi: 10. 10.11377 0137///.780887778833
[34] William E Schiesser. The numerical method of lines: integration of partial differential equations.
Elsevier, 2012.
[35] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In Proceedings of the 32nd International Conference on Neural Infor- mation Processing Systems, NIPS 18, page -552,, Red Hook, NY, USA, 2018. Curran Associates Inc.
12

---

[36] Yiheng Du, Nithin Chalapathi, and Aditi Krishnapriyan. Neural spectral methods: Self- supervised learning the spectral domain. The Twelfth International Conference on Learning Representations, 2024.
[37] Junho Choi, Taehyun Yun, Namjung Kim, and Youngjoon Hong. for parametric pdes without data reliance. Computer Methods in Applied Mechanics and Engineering, 420:1 2024.
Spectral operator learning
[38] Sam Bond- Taylo,, Adam Leach, Yang Long, and Chris G. Willcocks. Deep generative modelling: A comparative review of vaes, gans, normalizing flows, energy-based and autoregressive models. IEEE Transactions Pattem Analysis and Machine Intelligence, 44(11):7327-733 2022. doi: P01100//PAM..22111. 16668.
[39] Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:2210.02747, October 2022. doi:
10. 00448550/arxiv.221000774
[40] Shitong Luo and Wei Hu. Diffusion probabilistic models for 3D point cloud generation. In Proceedings of the IEEE/CVF Conference on Computer Vision and Puttern Recognition (CVPR), pages 2837-2845, June 2021.
[41] George Papamakarios, Eric Nalisnick. Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1-64, 2021.
[42] Arash Vahdat, Karsten Kreis, and Jan Kautz. Score-based generative modeling in latent space. In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 1 1128--11302 Curran Associates, Inc., 2021. URL
neuriss. http://sccccpppeffflleeeeeeee/ fiaer/ 200/211166cacb222222222266666113
[43] Arthur Bizzi. Matias Grynberg. Vitor Matias, Daniel Perazzo, João Paulo Lima, Luiz Velho, Nuno Gonçalves, João Pereira, Guilherme Schardong, and Tiago Novello. Flowing: Implicit neural flows for structure prruetureprreeerving morphing. arXiv preprint arXiv:2510.09557 2025.
[44] Yogesh Verma, Markus Heinonen, and Vikas Garg. Climode: Climate and weather forecasting with physics-ifformed neural odes. arXiv preprint arXiv:2404. /0024, 2024.
[45] John P Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.
[46] Claudio Canuto, M. Yousuff Hussaini, Alfio Quarteroni, and Thomas A. Zang. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. Springer, 2007. doi: 0000//9778------077----
[47] Lloyd N. Trefethen. Spectral Methods in MATLAB. Society for Industrial and Applied Mathe- matics, 2000. doi: 10.1 137/1. 717177997899999000
[48] Arthur Augusto Coutinho Bizzi. Flow-Structured Networks for Physics- Machine Learning. PhD thesis, Instituto de Matemática Pura Aplicada, Rio de Janeiro, Brazil, February 2025. URL hwt.h::tpp/.. h/ htttps///mpa htt/spppcootee///ppladd//////// /o9/outt/////t_____// _ese__
Arthur-Augusto- Cocttooooouuooooo-ooe Sitttinho-Bizzi.,. Bizzi. pddf
[49] Lawrence C. Evans. Partial Differential Equations, volume 19 of Graduate Studies Mathe- matics. American Mathematical Society, 2 edition, 2010.
[50] Arthur Bizzi, Lucas Nissenbaum, and João M Pereira. Neural conjugate flows: A physics- informed architecture with flow structure. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pages 15576- 15586, 2025.
[51] Marcelo Viana and José M Espinar. Differential equations: a dynamical systems approach to theory and practice, volume 212. American Mathematical Society, 2021.
[52] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary Vito,o, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in PyTorch. 2017. httpstts::/appi www.e/isaaaaiiicccciooemm
[53] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary Vito,o, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy,
13

---

Benoit Steiner. Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, high- performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché- Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, vol- ume 32. Curran Associates, Inc., 2019. URL
www.:ss..eu..pp..p ccppaper_ pdf
[54] Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Yee Whye Teh and Mike Titterington, editors, Proceedings of the Thirteenth International Conference Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Reseurch, pages 249-256, Chia Laguna Resort, Sardinia, Italy, 13-15 May mll
2010. PMLR. /tLL:ttts///p/oceeeiiii http://nws...l...ss/s/ss//////// Vttesss/v//gllootl11 glorot html
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In J. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neurul Information Processing Sys- lems, volume 20. Curran Associates, Inc., 2007. URL http::t::t:::tcttttttttt meers cc/ Paper pdf [56] Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, Jinkyoo Park, and Stefano Ermon. TorchDyn: implicit models and neural numerical methods in PyTorch. In Neural Information Processing Systems, Workshop on Physical Reasoning and Inductive Biases for the Real World, volume 2, 2021.
[55]
[57] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint ArXiv: 6980, 2014.
wrgg/bs////e 6980
[58] Sifan Wang, Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. An expert's guide to training physics-informe neural networks. arXiv preprint Stti::2308.08400 2023.
[59] Xinquan Huang and Tariq Alkhalifah. Pinnup: Robust neural network wavefield solutions using frequency upscaling and neuron splitting. Journal of Geophysical Research: Solid Earth, 127
2022. doi: https://doi.orr//00 10129//0//000000/2222011 10B0223700
[60] Shutong Qi and Costas D. Sarris. Hybrid physics-informed neural network for the wave equation with unconditionally stable time-stepping. IEEE Antennas and Wireless Propagation Letters,
2024. doi: 10. (099LAWP.2224.3355585
[61] Majid Rasht-Behesht, Christian Huber, Khemraj Shukla, and George Em Karniadakis. Physics- informed neural networks (PINNs) for wave propagation and full waveform inversions. Journal of Geophysical Research: Solid Earth, 127(5))e2221J0222) 1029/2021JB023.
2022. doi: https://doi.or//10
[62] Apostolos Parasyris, Lina Stankovic, and Vladimir Stankovic. Synthetic data generation for deep learning-based inversion for velocity model building. Remote Sensing, 15(11):2901, 2023.
[63] Aline Brougois, Marielle Bourget, Patriek Lailly, Michel Poulet, Patrice Ricarte, and Roelof Versteeg. Marmousi, model and data. In EAEG workshop-pracee ppractccl aspects of seismic data inversion, pages cp-108. European Association of Geoscientists & Engineers, 1990.
[64] Chengyuan Deng, Shihang Feng, Hanchen Wang, Xitong Zhang, Peng Jin, Yinan Feng, Qili Zeng. Yinpeng Chen, and Youzuo Lin. OpenFWI: Large-scale multi-structural benchmark datasets for seismic full waveform inversion, 2023. URL httppp///axxxv.orr/bbb... 02926
[65] Jesús Cuevas-Maraver, Panayotis G. Kevrekidis, and Floyd Williams. The sine-Gordon model and its applications. Springer, 2014. doi: 1111007/977-3331-00777
[66] Mahboub Baccouch. Numerical methods for the viscid and inviscid Burgers equations. In Computational Fluid Dynamics, chapter IntechOpen, Rijeka, 2024. doi: 10.5772iintcchppn. 1007351.
[67] Majid Khan. A novel solution technique for two dimensional Burger's equation. Alexandria Engineering Journal, 53(2):885-49. 2014. ISSN 1110-0168. doi: 10. aej6100ejjj.20144 .004.
[68] G. Batchelor. An introduction to fluid dynamics. Cambridge mathematical library. Cambridge University Press, reprint edition, 2005.
[69] Sebastian Acosta, Jesse Chan, Raven Johnson, and Benjamin Palacios. Pseudodifferential models for ultrasound waves with fractional attenuation. SIAM Journal on Applied Mathematics, 84(4):1609-1633 2024. doi: 137/2444411631
14

---

[70] Christiaan C. Stolk. A pseudodifferentil equation with damping for one way wave propagation inhomogeneous acoustic media. Wave Motion, 40(2):111-121, 2004, doi: 10.1016/j wavemoti. 1003.12.016.
[71] David E. Stewart. Numerical Analysis: A Graduate Course. CMSICAIMS Books in Mathemat- ics, 4. Springer, 2022.
[72] Nischay Rai and Sabyasachi Mondal. Spectral methods to solve nonlinear problems: A review.
Partial Differential Equations in Applied Mathematics, 4: 4:100043, 2021.
[73] Behnam Neyshabur, Zhiyuan Li, Srinadh Bhojanapalli, Yann LeCun, and Nathan Srebro.
Towards understanding the role of over- wee--aarmetrizationn in generalization of neural networks.
arXiv preprint arXiv: 1805. 12076, 2018
[74] Fabio Luporini, Mathias Louboutin, Michael Lange, Navjot Kukreja, Philipp Witte, Jan Hück- elheim, Charles Yount, Paul H. Kelly, Felix Herrmann, and Gerard J. Gorman. Ar- chitecture and performance of devito, a system for automated stencil computation. ACM Trans. Math. Sofnw., 46(1), apr 2020. ISSN 0098-3500. doi: 10.1145/3374916. URL https://dol org/ 1tlg//1..or 1145/3374916
[75] Felipe Linares and Gustavo Ponce. Introduction to nonlinear dispersive equations. Springer New York, 2015. doi:
A Causality of NeuSA
NeuSA is a "causal" architecture in the sense that it produces solutions which are evolutionary by nature, implying uniqueness and continuous dependence on initial conditions. It is generally impossible to prove the convergence of numerical method to the solution of a general PDE without imposing strong restrictions upon the spectrum of initial conditions and the function We may nevertheless prove that the solutions generated by our method have the properties associated with solutions to evolution problems. For the fnite-dimensionll system in eq. (4), these properties
are encapsulated in the properties of the associated flow operator T [0, 7] C X.xxxxXXX
C x Xxxxxxxx defined as:
t K'ún úo L. ((((T))dd
(14)
These operators have the following semigroup properties:
1. There exists an identity element: Φh "n ûo
2. The flow is an additive group action: f5' ûig- = 2b2 + un
Most importantly, these two properties translate immediately to the causal properties that we aim to prove: Property implies that initial conditions are enforced, while Property equivalent to uniqueness [51]. Theorem then follows as consequence of the fact that NeuSA solutions define flow.
Let Sb C (22) be the minee-dimeniional subspace spanned by our basis elements b. Consider now the decomposition operator P : Sb + C (o1xx--xexx an isomorphism mapping each element of S₁ into the tensor of coefficients of its expansion over the elements of b. We also define the reconstruction operator p Ca Xx,xx > Sp mapping the coefficient tensor onto the respective linear combination of the vectors bg- Inference for NeuSA consists of the following steps:
1. obtaining expansion coefficients ún Pug of the initial data in the basis b (possibly, projecting onto Sb);
2. acting on un by the flow of the vector field Fa using NODE:
3. reconstructing the continuous solution via u(t) Pfú(t).
These steps can be summarized as
p/ún úp ( F (((T)))
u (t) Pi
where
(15)
Now can formulate and prove the following
15

---

Theorem 2. For band-limited initial conditions uo E So and globalty-ipschitz neural vector fields Fo. the orbits created by NeuSA satisfy the initial conditions and are unique:
1. fulfillment of initial conditions: ug(0, x) x);
2. uniqueness: (0, -) (0, ub((,)) + ( VI [0,T]
Proof. For band-limited functions Sb. the decomposition P and reconstruction pi are bijections, and p is the inverse of P:
ptPu=
(16)
This fact immediately follows from uniqueness of expansion coefficients of any function E Sbb for the given basis b. Likewise, for a globally Lipschitz neural network Fo, there exists flow db associated with the NODE dů/dt (ŭ) (i.e., the flow determined by the vector field F Fo(ū)). Moreover, its orbits are unique. This follows from the classical existence and uniqueness of solutions for ordinary differential equations in Banach spaces [7] 51].
Property then follows from the uniqueness of the spectral decomposition combined with the flow properties of 4g:
Ug (0) = p Pu P =
(17)
Likewise, to prove property 2 we use the uniqueness of the spectral decomposition as well as the fact that NeuSA encodes the flow Do (under which the orbits do not intersect), represented as follows:
pt Pus
p PP 4/ Puo pt Puo us(12 + t1).
(18)
As mentioned above, in practice decomposition over the spectral basis is preceded by the projection of the initial data onto Sp. The basis can always be chosen to be large enough approximate any smooth function with the required accuracy. Likewise, it is generally reasonable to assume that common (finite-sized) feedforward networks are globally Lipschitz-continuous, since they essentially compositions of Lipschitz maps and Lipschitz nonlinear activations.
B Implementation details
In this section, we provide further details concerning our code implementation. The construction of the Fourier basis together with the extraction of its coefficients is explained in Subsection B.1 Lastly, Subsection dissccssss the numerical scheme adopted enhance dimension wise layers for multiple space dimensions, since they enable efficient modeling of complex interactions across spatial dimensions while avoiding the prohibitive costs of fully dense layers.
discusses
B.1 Spectral decomposition
To operate in the spectral domain, we represent the target functions using a spectral decomposition. The following subsections describe how the Fourier basis is constructed and how the corresponding coefficients are extracted.
B.1.1 Constructing the basis
We assume for this section that the domain Ω be decomposed direct product (e.g.. into direct product of 1D intervals), leading to the following product representation for the basis:
(x) = II by, (xi).
(19)
For example, the standard Fourier basis functions may be expressed the products of the following form (up to a normalization constant):
(x( x II exp(iwn; xi),
(20)
16

---

or, likewise, in terms of the respective sine/cosine functions, depending on the boundary conditions that have to be imposed.
In addition to enabling our initialization procedure, initializing b as the Fourier basis tends to result in dense representations for nonlinear and/or translation-invariin dynamics
72
While in the
context of numerical solution of PDEs this is often undesirable. in the context of neural networks this will become an advantage, denser connections between coefficients should lead (positive) overparametrization 73
Nevertheless, the Fourier basis has limitations. Foremost among them is its strictly non-local nature, which often results in difficulties when reproducing highly localized PDE solutions.
B.1.2 Extracting the coefficients ú and evaluating derivatives
In order to perform the spectral decomposition u + ú, we sample the function over N collocation N and represent it in terms of N basis terms. This may be expressed as:
points {x') ]1=0
u ((txx')
Σ (x(t))x(x'), i IE,,...,N,
(21)
where & is a dd-iimensional multi-index with N elements. This is a linear system that can be solved in straightforward manner, for example, given a grid structure, the Discrete Fourier Transform and its variations may be used. As it is usually done in the literature on spectral methods, we can simplify the basis representation by reshaping the domain via transformation X, leading to:
u(t, x²)
(x(t)bx((xx)))
(22)
Differentiation then takes place as described Section 2:
d ult, x²) Σ (t)(t) dx, bx(x(x')),
(23)
dx,
where the chain rule is used to calculate the right-most term. For a concrete example, 2D problem expressed in the Fourier basis over a square domain = | L, L) x [-L, L) may be mapped into the canonical domain |-x, T] x [-T,, T with the transformation xx) Tx/L, leading the following form for the derivatives:
u(t,x')
Σ
inTwk
II i=1
exp (iwx,
(24)
dx,
L
d dx3
ult, x(t,xx)= =->
Σ
Twky
(t)
exp II 1=1
exp (w.,,
(25)
In general, linear wranslation-inaarian differential operators may be obtained as polynomial on the frequencies w, leading to the following representation for the Fourier multiplier M used part of the initialization:
d² dx;dx;
F Finnca(u, (u, Vu, Vu,... := aou Σ
u
+
Σ 1.
aii
u
(26)
dxi
F lineur (i) = MO ú, with M₁ Go Σ O1s
iTTWks
Σ
a26
n mT Wk;
)
+ t..,, (27)
L?
ie. the k-indexed element of the multiplier M is given as a polynomial over the frequencies w. This multiplier may then be used for the initialization procedure described previously. For example, the Laplacian operator in space dimensions may be represented as:
π ?
π
? Kg
M
(28)
122
1,2
Analogous schemes also used for the sine and Fourier basis.
"This can be easily seen from the Fourier convolution theorem as appears when expressing nonlinear or non-trane www---nnaattoniivvaaeea invariant equations in the Fourier basis: pointwise products become convolutions involving coefficients that possibly far apart, resulting in dense connectivity matrices.
17

---

B.2 Dimension-wise layers for multiple space dimensions
Spectral methods inherently induce dense interactions between coefficients, even across distant modes. However, in multi-dimense dimensional settings, employing fully dense layers becomes computationally prohibitive. In particular, we can discuss the case of the 2D wave propagation problem presented in the experiments section: note that discretizing each dimension with just 100 spectral coefficients yields a total of 104 basis elements. single dense linear layer operating on this space would then require (110)) parameters. As a result, when such layers are applied repeatedly (as would be done by a Neural ODE using for vector field), the computational cost becomes excessive.
Similar challenges in computer vision have led to the development of convolutional neural networks (CNNs), which employ finite-support kernels to create localized linear connections among neigh- boring nodes. However, this approach is ill-suited to our setting, First, CNNs are built to embed the translation invariance inherent to image recognition tasks, which is a property that does not generally apply to our problems. Moreover, their local receptive fields restrict long-range interactions, which are essential in spectral representations.
Instead, it is necessary to construct a layer that performs the dense and global-reaching connections associated with spectral methods while remaining parameter-light. We have opted for low-rank, dimension-wise linear layers, as shown in Figure 8
These networks consist essentially of Hadamard (element-wise) product followed by alternating dense linear transformations applied row-wise and column-wise. This structure can be efficiently implemented by combining tensor transpositions with batched matrix multiplication operations available in PyTorch and similar deep-learning frameworks (see Algorithm [T. For instance, given 2D coefficient matrices útt) I m and layer parameters A Rmxn B Raxn, and C E mm this yields:
1. Element-wise scaling (Hadamard product):
ú + ú A A
where each element ü, is multiplied by the corresponding Atj.
2. Row-wise linear transformation:
T1
ri
T2
BBB
:
[mm
I'm B
where each row vector r; is multiplied by the matrix B.
3. Column-wise linear transformation:
[c
62
c,] + ICes Coz
Cenl
where each column vector ej is multiplied by the matrix C.
The matrix ú can be visualized
Ú1.2
Ú1,
02
02.2
12
T2
C2
Gn]
üm
üm.2
úm...
Lmm_
Each of these layers applied a 2D spatial domain of dimensions m, n requires 0(m)) parameters and operations, in contrast to the O(m?n' (mm2n)) complexity of a naive fully dense layer. This approach allows us to leverage the "densification" inherent to spectral methods, enhancing deep learning performance, while reducing computational and memory demands. In practice, we have found these layers to be effective drop-in replacements for dense linear layers when handling large, multi- dimensional inputs.
18

---

Algorithm Implementation of dimension-wise layers for two dimensional input
input ()
ú E R m
i < ú A
> Hadamard (pointwise) product, A E IR
< linear [ú, B)
Batched matrix multiplication, E Raxx
i < i²
Transpose > Batched matrix multiplication, C "
u - linear(tn (ü, C)
û < ú²
Transpose Nonlinear Activation
- act(ú)
As illustrated by comparing Figures and g these layers stand in contrast to convolutional layers, yet serve complementary purpose. While convolutional layers impose sparsity through local connectivity and translation invariance, dimension-wise linear layers aim to preserve sparsity while retaining long-ranging connections. They can also be interpreted as low-rank tensor applications or Kronecker products.
Figure 7: Schematic for convolutional linear layers.
Hadamard
Row wise
Column wise
Figure 8: Schematic for dimension- wisse linear transformations.
C Experimental details
In this section, we provide further information concerning hardwares and licensing. We also discuss additional details on the numerical experiments presented in Section 3 of the paper.
Hardware and resources. All experiments were executed on identical machines, containing an Nvidia RTX 4090 GPU with 24GB VRAM, Intel i9- 19-00KK processor, and 128GB RAM.
Licenses. The code is implemented in Python using the libraries PyTorch (version 2.1.0) and torchdyn (version 1.0.6), distributed under the BSD 3-Clause and Apache licenses, respectively.
General setup for numerical experiments. All experiments are evaluated against ground-truth solution obtained using a state-of-thed numerical method. For each model and experiment, we report the emor metrics rMAE (relative Lı error) and rMSE (relative L2 error) between the predicted solution Upred and the ground-truth solution uGT, computed over N evaluation pairs (4,, x;), as follows:
Σ [upred (li, xs) - uCT(li, x:)
iMAE
(29)
ΣN [UGT(tixx))
ΣN (upred (li, xi) u (lis x;))²
rMSE
(30)
Σ (uer(l,, x²))²
19

---

we can rewrite the problem as a system of In the particular case of the Burgers' equation, where the output is a 2D vector (u, v), the rMAE and rMSE are computed separately for each component, and their average reported the final error metric.
Due memory constraints, it is feasible feed the entire spatial grid the MIP--based architectures for 2D equations. Instead, we adopted random uniform sampling at every step. For PINN, QRes, and FLS models, we sample 10, 000 points for the PDE residual, 1000 points for the initial condition and 500 points for the boundary condition. In contrast, PINNsFormer generates temporal sequence of five collocation points for each sample, which are then processed through its pecoder-decoder architecture. To ensure comparable memory usage and training time with the other baseline methods, we train PINNsFormer on the Burgers' and wave equations using 2, 000 points for the PDE loss, 200 for the initial condition, and 100 for the boundary condition. This corresponds to one-fifth of the amount sampled for other baseline methods.
C.1 Wave equation
Ou O
For the wave equation (TT. by introducing first-order equations (with respect to time) given by
0t u1 V,
(31)
Ot 1 = (x)Au.
The vector field to be learned by NeuSA is defined as:
6
Fo(ū) = M Mon+
(32)
where (ù, i), the weight and the entries of the matrix M as defined in (88) Using simplified version of Fg that takes only à as imput accelerates learning, as 0 does not influence the equation, and also enables more compact model. The velocity fields may be visualized below, in Figure
-2.0
1.5
15-
1.4
-1.0
0.5
0.0
1.2
1.0
11
1.5
-25
2.0
-1
1.0
-20
-2.
-1.
-1.
nb
1.6 2.0
(a)
Figure 9; The 2D three-layer and Marmousi mediums.
To train the NeuSA model, we simulate an infinite domain by extending the original spatial region from 2| to [-4,4] The extensions of the original domain are designed to ensure that, within simulation's 's temporal window, any reflected waves do not re-enter the original region of interest. This effectively prevents boundary artifacts from interfering with the solution. We use 201 d basis elements and integrate over 201 time steps, creating grid with spatial step of Ar Ay = 0.04 and temporal step of At 0.01. For the spectral decomposition, adopt cosine basis.
20

---

Figure 10: Solutions for the 3D wave experiment.
For the baseline models, we set the initial condition weight to 10 (on both Dirichlet boundary condition and first-order initial condition, see Eq. (11) in the main text). This weighting led to consistently improved accuracy across experiments.
The ground-trui solution is calculated using standard finite central difference scheme for the second derivatives with order 8 in the space and order 2 in time (74), with spatial step of 0.01 and time step of 0.001. The domain was also extended simulate infinite domain.
C.2 1D sine-Gordon equation
du
Let us define again in the sine-Gordon equation 12) which allows us rewrite it as an evolutionary system of the form
a
Öt u C,
(33)
0
V
Au 10 (u).
Øt
From this and Eq. (12), the vector field to be learned from NeuSA can be defined as:
(ū) ( (
6 ( cF(())
(34)
where the weight is set to 0.1, and the entries of the matrix M are taken as defined in 28
Owing its architectural design, NeuSA automatically satisfies the initial condition, while the use of a sine basis ensures compliance with the boundary condition. For the baseline models, we adopted a weight of 1 to the initial condition, this choice yielded the most accurate results in our preliminary experiments.
In contrast to the experiments involving 2D equations, training the baselines for the sine-Gordon equation can be performed using the full spatial grid at each step. Since PINNsFormer produces training points per grid location, is trained
coarser regular grid of size 101 101, whereas all other models are trained on a finer 201 x 201 grid. This choice of grid size presented a reasonable balance of training time, accuracy, and memory usage.
The ground-truth solution was computed through pseudo-spectral method combined with 4th-order Runge-Kutt integrator. We used time step of At 0.001 and spatial resolution of Ar 0.04.
21

---

C.3 2D Burgers' equation
Two-dimensional Burgers' equation for the vector function (u(t, z, v), v(t, 22,, y)) can be written in the form of system
0 vāu
du
du
u =
u
V
01
8x
0y
(35)
0
0v Ov
vAv
01
dr
Dy
From this and (13), the vector field for NeuSA can be directly defined by
PM o u ( Fo = vM 1 (0t
(36)
where M is given in (28). and the weight is set to 0.1. Also, (ú, 0) and (ū, 0) are the neural networks whose parameters we optimize during training, so the loss function is composed by the residues of each of them. The equal subscript 0 on both is slight abuse of the notation, since the parameters of these networks are not the same.
NeuSA inherently satisfies the initial conditions by construction (Theorem 1), and the use of Fourier basis the Burgers' equation ensures that periodic boundary conditions are also met automatically. For the baseline models, we set the initial and boundary conditions weights to èrc = 10a and λBC = 1, as these values produce the best performance in preliminary experiments.
The results presented in Table for the Burgers' equation correspond to model training and evaluation over the time interval [0, 1]. We compare them and the extrapolation experiment with ground-truth obtained through pseudo-spectral method integrated with 4th-order Runge-Kutta for the time interval [0, 2]. The temporal and spatial discretizations are to At 0.001 and Ar Ay 0.02, respectively.
C.4 Quantitative results
The 2D experiments (Burgers' and wave equations) were run with different seeds each, while the ID sine-Gor Gordon experiment was run with different seeds. The relative L₁ and L2 metrics and training times (in seconds) presented in the paper are an average over all random seeds. Tables and 6) we list the mean and standard deviation of these values for the main experiments reported in the paper.
Table 2: 2D Layers: mean and standard deviation of rMAE, rMSE and runtime (TT).
Model
rMAE
rMSE
TT
mean
std
mean
std
mean
std
PINN
7.66 10
1.27 x 10-1
5.45
10-1
8.31
10
566
7.2
QRes
1.54 x 10
3.45 x 10
1.15 x 10
3.10 x 10
750
2.4
FLS
8.87 x 10-1
1.63 x 10-1
5.90
10-1
9.96
10
577
5.5
PINNsFormer
1.56
2.78 x 10-1
1.07
1.41
10-1
1,484
3.4
NeuSA (ours)
1.02 x 10-1
3.42 x 10
7.49 x 10
2.24
10
530
1.4
Table 3: Marmousi: mean and standard deviation of rMAF, rMSE and runtime (TT).
Model
rMAE
rMSE
TT
mean
std
mean
mean
std
PINN
1.03
1.48 x 10
6.98 x 10-1
9.48 x 10 2
635
29.7
QRes
5.57 10
6.03
x 10
4.12 x 10-1
3.06 x 10 2
718
15.8
FLS
9.99
10
1.50
10
6.84
10-1
9.02
10
648
23.9
NeuSA (ours)
2.20 x 10-1
4.43
10
1.71 x 10-1
3.45 x 10
573
2.28
22

---

Table 4: 3D Layers: mean and standard deviation rMAE, iMSE and runtime (TT).
Model
rMAE
rMSE
TT
mean
std
mean
std
mean
std
PINN
1.86 x 10
6.07 x 10
7.32 x 10 2
1.75 x 10 2
5990
37.3
QRes
4.23 10
7.40
10
2.12
10
2.60 x 10
8775
30.2
FLS
1.70 x 10-1
4.46
10-2
6.95
x 10-2
1.32 x 10
6179
41.0
NeuSA (ours)
1.18 x 10
2.50 x 10-3
8.40
x 10-3
1.80 x 10-3
702
7.3
Table 5: ID Sine-Gordon. mean and standard deviation of rMAE, rMSE and runtime (TT).
Model
rMAE
rMSE
TT
mean
std
mean
std
mean
std
PINN
1.70 10-1
3.31
10-2
1.39
10-1
4.70
10-3
976
41.4
QRes
2.58 x 10-2
4.50
10-3
1.99
10-2
2.80 x 10-3
1,315
49.0
FLS
1.43 x 10-1
7.20
10-3
1.35 x 10-1
2.01 x 10-2
1,015
68.3
PINNsFormer
7.85 10-1
4.95
10
6.81
10
3.31
10
3,333
186.6
NeuSA (ours)
1.23 x 10-3
1.49
10-5
9.16
10-4
5.98 x 10-6
215
1.3
Table 6: 2D Burgers: mean and standard deviation of rMAE, rMSE and runtime (TT).
Model
rMAE
rMSE
TT
mean
std
mean
std
mean
PINN
1.65 x 10
8.91
10
2.21 x 10
9.30 x 10-2
871
2.9
QRes
2.32 10
1.90
10-3
7.27
10
1.60
10-3
1,5
2.3
FLS
1.47 x 10-1
7.68
10
2.02
10-1
8.08 x 10-2
885
0.3
PINNsFormer
1.06
1.91
10
1.05
1.71 x 10
2,294
108.1
NeuSA (ours)
7,66 10-2
1.96
10
1.15 x 10
1.71
10
62
0.2
D Additional experiments
D.1 Number of frequencies and Training time
Due to its built-in integration process, NeuSA's convergence depends on guaranteeing the numerical stability of the Neural ODE. particular, the more basis elements are used, generally the smaller the time-steps for the integration should be. This introduces a trade-off between spatial resolution and speed. We exemplify this behavior through the following additional experiment: We consider the Burgers Equation, with initial conditions similar to those in the main paper, and evaluate how long takes to train NeuSA with a varying number frequencies. The results can be seen Table
Table 7: Number of frequencies vs rMSE and training time. More elements are more accurate, but slower.
Number of frequencies
rMSE
Training time (s)
61
7.6e-2
109
81
5.6e-2
146
101
4.3e-2
183
121
3.3e-2
223
141
2.6e-2
264
23

---

D.2 Inference time
As discussed, NeuSA's inference time is dominated by the process of integrating the Neural ODE. In essence, inference and training of Neural-Spectral Architectures is relatively costly, due to the extreme effective depth of NODEs; on the other hand, one pass generates the solution at all space and time points. Nevertheless, the much faster convergence rate in terms of number of epochs generally offsets this effect.
We compare the inference time for a subset of our experiments in Table 8 where we evaluate all models over uniform space-time grid. As can be seen from the table, NeuSA's inference time remains more or less constant across the experiments. In contrast, the baseline architectures process every point in the spatio-temporal grid simultaneously, resulting in substantially slower inference when iterating over large amounts of points, likely as a result of memory swapping. Throughout our experiments, we attempt to mitigate this effect for the reference architectures by randomly sampling the domain in the experiments with the 2D Wave/2D Burgers equations.
Table 8: Time in seconds required by each architecture to compute the solution at the grid-like like collocation points used in the experiments.
Architecture
Wave 2D (101 x 101 x201)
Burgers 2D
Sine-Gordon (201 x201)
MLP
0.04
0.28
0.0005
QRes
0.11
0.34
0.002
FLS
0.04
0.28
0.0006
PINNsFormer
1.29
2.57
0.02
NeuSA
0.12
0.11
0.10
E
Limitations
As mentioned in Section 2.1, NeuSA requires the initialization of the linear model, which may not be obvious. For Burgers' equation, for example, we have obtained the best results initializing NeuSA near solution the associated heat equation. As it is well known, however, some linear PDEs may turn out to be stiff when expressed in the spectral domain. This stiffness may introduce instability into the integration of the Neural ODE, For example, when we apply spectral decompositions to the Korteweg-de Vries (KdV) equation 775 the resulting system of ordinary differential equations becomes essentially stiff, especially when wide range of frequencies is considered. This will lead to excessively large eigenvalues of the matrix used in the implementation of the associated solver, and may result in a numerically unstable solution.
This numerical issue, however, does not affect a large number of equations important for various practical upplications. Consider, for instance, the canonical example of the wave equation:
Au
(37)
Transforming it into the first-order hyperbolic system as shown above and applying Fourier transform arrive at the following ODE system
d ④ dt
0 D
! and D
A
where A
(38)
i
The block matrix A has eigenvalues X x= iik for A: == 1,..,, N. According to the standard stability criterion for the Runge-Kurte Kutta methods (based on the stability polynomials and eigenvalues of the equation matrix), time step of h O (+) is required for the latter system, which is significantly larger than the time step required for the ODE system associated with the KdV equation.
This indicates that the Runge- Kutta method allows for larger time step when applied to the wave equation compared to the KdV equation, thereby improving the practicality and efficiency of our method for broad class of problems, including those involving wave phenomena.
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The paper proposes new of architectures, called Neuro-Spectie Spectrlll Archi- tectures (NeuSA), which are clearly described and delineated Section 2. The properties claimed about this architecture experimental setup through which we verify these claims, while the results are finally concluded in Section 4, and consistent with abstract and introduction.
explained in the same section. Section 3 presents the
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize other settings.
is fine to include aspirational goals as motivation as long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The limitations of the work performed are mentioned in Section 4. Details about our limitations as well as some experimental results showing these can be found in Appendix D.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model vell-specificattin, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
25

---

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

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Question:
Answer: [Yes]
Justification:
Code will be provided via GitHub for the camera-ready version of the paper.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ //nips fliaps.cccpublice wwtppbbiccgguueeeeeoooooocppcped
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The hyperparameters and optimizers have been clearly indicated. All in- formation for running the experiments can be found in Section and Appendices and C.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[Yes]
Justification: The experiments were reported with an average over multiple seeds, with standard deviations reported in Appendix C.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emoo bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
27

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative eπτor rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
The computer resources used for the experiments have been clearly written.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics htttp:://eeurjs: cc/ppublcc/thiiccciiidiee....
Answer: [Yes]
Justification:
The research does not in any way conflict with NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer:
[Yes]
Justification: As pointed out in the introduction and the experiments, NeuSA well suited for many physical applications in areas such as biology, cosmology and physics. This may lead positive impact in any many applications related to physical modeling.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
28

---

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

---

If assets are released, the license, copyright information, and terms of use in the thcode. has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssi
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA].
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA].
Justification: The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
30

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA].
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM for what should or should not be described.
31


Reviews:
[{"review_id": "fz8zpqW2sV", "reviewer": "Reviewer_3zAY", "summary": "This study introduces NeuSA (Neuro-Spectral Architectures), which is an adaptation of the Neural Ordinary Differential Equation (Neural ODE) architecture to the Physics Informed Neural Network (PINN) problem - that is, approximating the solution of a given PDE by using the PDE residual in the loss function. Inspired by the spectral methods for PDEs, NeuSA models the solution as a linear combination of spatial basis functions and uses a Neural ODE to learn the time dependence of the basis coefficients. Furthermore, while the basis functions are initialized to be the Fourier basis, these are also learned during training by the use of learnable orthogonal transforms. As the basis function expansion can be set up to satisfy both the initial and boundary conditions, the loss function only contains the PDE residual term. Experiments on challenging low dimensional systems show the strong performance of NeuSA compared to the baselines considered.", "questions": "The first part of my questions are related to the points raised in the Weaknesses section above.\n\n1. How does the algorithm performance and the computational overhead change with the number of basis elements? Is this a hyperparameter to be sweeped over, or do the authors have a strategy / heuristic to determine the appropriate number? \n\n2. Related to the previous question, the initial conditions for the PDE problems in the experiments are either well described with finite Fourier basis elements (2D Burgers': product of sines/cosines; 2D wave, 1D sine-Gordon: Gaussian, whose spectrum decays as $\\exp(-k^2)$ with increasing wavenumber $k$.). When the initial condition becomes more difficult to describe with a truncated Fourier series (such as a square wave), does this method still retain its performance?\n\n3. I am curious about the per-epoch training time for each of the models considered in the paper. Is my reasoning correct, or have I made a mistake?\n\nSome more minor questions are as follows:\n\n4. Neural ODEs tend to struggle especially for long time series data, as the gradient pathologies become increasingly likely to occur (due to all the recurrent computation that stacks up while propagating the ODE solver forward in time) and the cost of the algorithm also scales linearly with time. In the context of NeuSA, this situation would correspond to solving a PDE over a longer temporal domain. Is this observed in practice? \n\n5. The algorithm's initialization scheme relies on linearizing the PDE and doing some algebra to derive the Fourier multipliers. Was this done manually or automatized, for example, by symbolic calculations?", "limitations": "Yes", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "**Strengths**\n\n1. Strong performance on challenging problems\n\nThe problems considered in the paper such as wave propagation in heterogeneous media or 2D Burger's equation with shock development are all difficult problems that conventional PINNs tend to struggle on. However, the proposed approach returns accurate solutions for these problems with much less computational cost than its competitors.\n\n2. Stable training dynamics due to improved architecture design\n\nDue to the model formalism, NeuSA respects causality and does require multiple competing terms in its loss functions. This leads to faster and more stable training, as evidenced in the right panel of Figure 1.\n\n\n3. Clear architectural motivations\n\nThe use of a Neural ODE-like architecture for the PINN problem is clearly motivated from the viewpoint of spectral methods for PDEs.\n\n\n**Weaknesses**\n\n1. Choosing the number of basis function elements\n\nThe number of Fourier basis elements used for training seems to be a key hyperparameter for this method, as it affects multiple aspects of training, including: 1) the maximum spatial frequency that the model output can have, 2) the computational cost of the model ($n$ basis elements for each of the $d$ spatial dimensions will scale roughly as O($n^d$), even with the dimensionwise layers used in this study), 3) how well the initial condition is described by the model. As such, an experiment to gauge the model performance as a function of number of basis elements used would provide more insights about the proposed method.\n\n2. Confusion about the speed of the algorithm\n\nThe authors claim that despite the slowness of the Neural ODE architecture, NeuSA can achieve comparable or faster training times than competing non-Neural ODE approaches, with one of the reasons being NeuSA being able to compute the loss with a single forward pass whereas conventional approaches requiring multiple passes. While the shorter training times reported in Table 1 and Table 2,3,4 (in the Supplementary) seem to support this claim, it needs to be noted that the total training steps for NeuSA were 10 (2D Wave, 1D Sine-Gordon) or 100 times lower (2D-Burgers) than its counterparts. While impressive, this does mean that the per-step training time of NeuSA is often much higher than the other approaches (for 2D Wave, x3.57-x9.40 times slower), and the above claim does not seem to hold.\n\nMinor\n\n3. While recently proposed, a method for solving PDE that draws from the spectral method does exist in the literature [1]. While that approach and the current paper differs quite a bit, it would be great if the authors could incorporate the reference in their discussions, especially as the reference considers non-Fourier basis functions as well.\n\n- [1]: J. Choi et al. Spectral operator learning for parametric PDEs without data reliance. Comput. Methods Appl. Mech. Engrg. 420, 15 (2024).", "was_revised": true, "final_justification": "This is a solid paper on developing a Neural ODE-based PINN architecture that shows strong performance on challenging PDE problems. The concerns I had prior to the discussion phase were minor, and even these were well-addressed by the authors with multiple extra calculations performed. As such, I will keep my original positive score.", "rebuttal": "We thank the reviewer for their endorsement and detailed analysis. We are encouraged that they find our architecture to be well motivated and interesting. We appreciate that the reviewer highlights the model's strong experimental performance and notices its substantially faster convergence speed. We respond in detail to the weaknesses and questions raised by the reviewer below. Weaknesses 1 and 2 are addressed in questions.\n\n## Weaknesses\n\n> 3. **While recently proposed, a method for solving PDE that draws from the spectral method does exist in the literature. While that approach and the current paper differ quite a bit, it would be great if the authors could incorporate the reference in their discussions [...]**\n\nWe appreciate the additional reference and will include it in the final paper. As mentioned, the approaches are quite different: The authors of SCLON are interested in operator learning for solving parametric partial differential equations (PDEs), while we address a fixed initial value problem, more in line with PINNs literature. \n\nAlthough they also leverage expansions in orthogonal functions, the way that time propagation is dealt with is also quite different: SCLON divides the time domain into multiple discrete subintervals and trains a separate network for each; our method, in turn, extensively uses continuous-time Neural ODEs. Nevertheless, some aspects — such as Legendre spectral methods and the sequential approach for time-dependent PDEs — could be incorporated into future extensions of the NeuSA framework.\n\n## Questions\n\n> 1. **How does the algorithm performance and the computational overhead change with the number of basis elements? Is this a hyperparameter to be sweeped over, or do the authors have a strategy / heuristic to determine the appropriate number?**\n\nOur initial criterion was to use as many basis elements as possible for a given spatial sampling of the domain. To promote a wider discussion on the role of basis functions, we have included a small section in the supplementary material with experiments exploring this question further. \n\nWe consider the 1D Burgers equation with a variable number of basis elements. To ensure numerical stability, we increase the number of time steps proportionally to the number of basis elements.  The table below shows how the relative root mean-square error and the training time are affected by the number of frequencies used.\n\n| Number of frequencies | rMSE | Training time (s) |\n| -------- | ------- | ------- |\n| 21 | 0.181 | 35 |\n| 41 | 0.109 | 72 |\n| 61 | 0.076 | 109 |\n| 81 | 0.056 | 146 |\n| 101 | 0.043 | 183 |\n| 121 | 0.033 | 223 |\n| 141 | 0.026 | 264 |\n| 161 | 0.021 | 299 |\n\nAs can be seen, the number of frequencies has a substantial influence over both the error and training time, but the accuracy benefits have diminishing returns. As a general heuristic for classical spectral methods, we would expect small improvements after the basis elements chosen cover the frequency band induced by the solution; it is still unclear if or how this intuition carries over to the neural case. We will include this experiment in an appendix.\n\n> 2. **Related to the previous question, the initial conditions for the PDE problems in the experiments are either well described with finite Fourier basis elements [...] When the initial condition becomes more difficult to describe with a truncated Fourier series (such as a square wave), does this method still retain its performance?**\n\nThe effects of high-frequency components on performance may be seen in the results of Experiment number 3, with the 2D Burgers' equation. The solution eventually converges to a sharp, discontinuous shock profile with high frequency content, which is widely known to be a difficult scenario for spectral and pseudo-spectral methods due to the non-localized behavior of the Fourier basis elements. Nevertheless, NeuSA ranks second among the listed architectures, while vastly outperforming the others in terms of speed.\n\nIn order to investigate NeuSA's performance on initial conditions less easily described by the Fourier basis, we have performed an additional experiment on the Sine-Gordon equation, with a triangular initial condition, the quantitative results of which are shown in the table below:\n\n| Model | rMSE | Training time (s) |\n| -------- | ------- | ------- |\n| PINN | 4.47E-02 | 420 |\n| QRes | 1.27E-02 | 498 |\n| FLS | 4.43E-02 | 433 |\n| NeuSA | 2.02E-02 | 292 |\n\nWe observe that NeuSA has a slightly worse performance than QRes, albeit with a smaller training time. This can be attributed to the slower decay of the triangle wave's Fourier coefficients.\n\n\n> 3. **I am curious about the per-epoch training time for each of the models considered in the paper. Is my reasoning correct, or have I made a mistake?**\n\nWe discuss this briefly in the paragraph starting at line 273:\n\n*\"Although built on the inherently slower NODE framework, NeuSA trains significantly faster than MLP-based PINNs\"*\n\nIn essence, inference and training of Neural-Spectral Architectures is indeed slower, due to the extreme effective depth of NODEs.  On the other hand, the much faster convergence rate in terms of number of epochs generally offsets this effect. \n\nFor the paper's final version we will include an additional snippet expanding on this point, as well as a table with the **average inference time (in seconds)** for each model. For the table, we include the time in seconds needed for each model to calculate the solution at the grid-like collocation points used in the experiments:\n\n| Architecture | Wave 2D (101x101x201) | Burgers 2D (201x201x101) | Sine-Gordon (201x201) |\n| -------- | ------- | ------- | ------- |\n| NeuSA | 0.12 | 0.11| 0.10 |\n| MLP | 0.04 | 0.28 | 0.0005 |\n| QRes | 0.11 | 0.34 | 0.002 |\n| FLS | 0.04 | 0.28 | 0.0006 |\n| PINNsFormer | 1.29 | 2.57 | 0.02 |\n\nAs can be seen from the table, NeuSA's inference time remains more or less constant across the experiments, as it is dominated by the numerical integration process. In contrast, the baseline architectures process every point in the spatio-temporal grid simultaneously, resulting in slower inference when iterating over large amounts of points. Throughout our experiments, we mitigate this effect for the reference architectures by randomly sampling the domain in the experiments with the 2D Wave/2D Burgers equations.\n\n> 4. **Neural ODEs tend to struggle especially for long time series data, as the gradient pathologies become increasingly likely to occur (due to all the recurrent computation that stacks up while propagating the ODE solver forward in time) and the cost of the algorithm also scales linearly with time. In the context of NeuSA, this situation would correspond to solving a PDE over a longer temporal domain. Is this observed in practice?**\n\nYes, we observed this in some experiments with very long time intervals. There are several ways to address this issue, but one approach that has been successful in our trials is to train different models over different time intervals. This increases the number of parameters and training time, but makes the model more adaptable to the different behaviors of the solution over different time intervals. Still, conventional architectures would generally fail under the same conditions, the longer time horizons inducing failure modes due to their lack of causality.\n\n> 5. **The algorithm's initialization scheme relies on linearizing the PDE and doing some algebra to derive the Fourier multipliers. Was this done manually or automatized, for example, by symbolic calculations?**\n\nThis was done manually. We will state this more clearly in the final version."}, {"review_id": "863mL9Sro2", "reviewer": "Reviewer_gAMj", "summary": "The paper introduces NeuSA (Neurospectral Architectures) as a novel class of PINN architecture. This class satisfies by construction the boundary/initial conditions of the PDE, which alleviates one of the main training difficulty of PINNs, given by the imbalance in the loss components. Moreover, time dependency is resolved through continuous time integration based on a NeurODE formulation. The model also satisfies causality and uniqueness of the solution, which are two very interesting properties for a Physics-Informed model. The proposed NeuSA is validated on some benchmark PDEs and compared with state-of-the-art architectures.", "questions": "* I find it surprising that QRes outperforms NeuSA in the relative L2 Error over time (except far from the extrapolation region) in Figure 5, do you have some intuition on why this happens?\n* I was wondering whether the authors tested also other training protocols (as quasi newton methods) which might enhance the performance of the Model in the training region. \n* Related to the question above, I would be curious to know (if the authors have tested it) what would be the extrapolation performance of the model when a NeuSA is trained to very good accuracy (say, around 10^-6 error) in the training region. Would the extrapolation error increase or stay limited?", "limitations": "The authors have addressed the limitations throughout the paper.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "The paper is in a very good format. In particular:\n* The model is well presented, and the architectural choices are well motivated. Overall, the paper is well written.\n* I find the combination of the PINN framework with the concept of NeurODE for the time integration very elegant.\n* The model tackles by construction several long-standing issues of the PINN architecture.\n* The achieved performance in training time is surprising, considering the time integration chosen.\n\nOn the other hand, there are some minor issues that concerns me:\n* The PDEs considered for testing seem very suitable for the chosen architecture, due to the choice of a Fourier basis.\n* The model only considers simple geometries, while I think the real challenge in imposing boundary conditions typically lies in more complex ones", "was_revised": true, "final_justification": "The discussion with the authors confirmed my score.", "rebuttal": "We thank the reviewer for their insightful observations and endorsement. We appreciated that they find our architecture to be well motivated and interesting, and noted the substantially faster convergence speed of our method. We are also glad that they acknowledged the architecture's unique features, including uniqueness of trajectories and causal structure. We address in the weaknesses and questions pointed by the reviewer below.\n\n## Weaknesses\n\n> 1. **The PDEs considered for testing seem very suitable for the chosen architecture, due to the choice of a Fourier basis.**   \n    \nThough we understand the reviewer's concern, we would argue that the PDEs chosen for the experiments are not particularly tied to the strengths of the Fourier basis. In particular, both the wave equation in inhomogeneous media and Burgers' equation involve highly localized phenomena, shocks and high-frequency components, which Fourier series struggle to represent.  \nTo illustrate this point, we will add to the final version an additional experiment on the 2-D wave equation, considering the Marmousi model [1]. This is a highly irregular acoustic model that induces highly complex reflections. This problem is a canonical benchmark in acoustic geophysics and is extremely challenging for even specialized PINNs [2].\n\nWe tested NeuSA on the 2D wave equation using a cropped version of the Marmousi model. By doubling the frequency content and training for 3000 epochs, we achieved an RMSE of 9.05e-2 with a training time of 923 seconds. This result is comparable to the framework tackled by [2], whose architecture is specialized to simulating seismic wavefields for nonsmooth complex media. This model has a much larger set of highly optimized hyperparameters and was trained using time-domain decomposition over 45000 epochs with 80000 collocation points, and requiring hours of training current setup. We therefore conclude that NeuSA has a strong performance in this scenario when compared to most of the literature. We will add these results to the final version of our paper.\n\n> 2. **The model only considers simple geometries, while I think the real challenge in imposing boundary conditions typically lies in more complex ones.**  \n\nWe thank the reviewer for their observation. We quickly mention this in the paragraph starting on line 283: \n \n*\"all experiments take place on simple spatial domains [...] more complex geometries would require the implementation of more generic bases,\"*  \n\nIn other words, NeuSA deals most easily with simple domain shapes, whereby we may leverage the basis's properties. This may be extended in a straightforward manner by using differentiable implementations of other bases (such as cylindrical and spherical harmonics).  \nOn the other hand, NeuSA remains a meshless PINN, so it is still possible to enforce different boundary conditions and domain constraints through carefully crafted losses.    \nWe also stress that many developments in PINNs, including the benchmark paper PINNacle [4] and the reference papers for PINNs, FLS, QRes, and FLS, deal almost exclusively with domains with simple boundary shapes.  \n\n## Questions\n\n> 1. **I find it surprising that QRes outperforms NeuSA in the relative L2 Error over time (except far from the extrapolation region) in Figure 5, do you have some intuition on why this happens?**  \n\nOur intuition is that this is due to the solution of Burger's equation developing shocks over time, which Fourier basis struggle to represent well. Due to its reliance on spectral bases, NeuSA cannot truly \"focus\" on any given part of the domain, unlike other architectures. In other words, the Fourier basis will tend to induce global components instead of local components of the solution. We have verified that the vast majority of the error is concentrated in the close vicinity of the shocks; NeuSA seems to approximate the solution quite well outside of this interval, as can be seen from the experiment results. \nOn the other hand, we note that QRes is a highly nonlinear and flexible network, and is able to better represent solutions with shocks. \n\n> 2. **I was wondering whether the authors tested also other training protocols (as quasi newton methods) which might enhance the performance of the Model in the training region.**  \n\nTo address this question, we have run new experiments where we use the LBFGS optimizer [5], instead of ADAM, and compared training time and rMSE. The results are in the table below\n\n| Equation | LBFGS rMSE  | ADAM rMSE | LBFGS training time (s) | ADAM training time (s) |  \n| -------- | ------- |  ------- |  ------- |  ------- |\n| Burgers 2d | 0.142 | 0.11 | 576 | 62 |   \n| Wave 2d | 0.0502 | 0.07 | 625 | 530 |  \n| Sine-Gordon | 0.00105 | 0.0009 | 831 | 215 |  \n\n As it can be seen in the table training with LBFGS often led to much slower training in exchange to little to no improvements in accuracy. We will add this table and discussion to our appendices in the final version.  \n> 3. **Related to the question above, I would be curious to know (if the authors have tested it) what would be the extrapolation performance of the model when a NeuSA is trained to very good accuracy (say, around $10^-6$ error) in the training region. Would the extrapolation error increase or stay limited?**\n\nTo address this question, we will add an additional experiment focusing on the performance of NeuSA for extrapolation. We trained the Sine-Gordon equation over the time interval $[0, 1.5]$ using 100 and 1000 time steps, and obtained PDE losses in this time interval of 6e-4 and 8e-6, respectively. The RMSE over the interval $[0, t_f]$ is shown in the table below:\n\n| Time Steps | Loss | $t_f = 1.5$ | $t_f = 1.6$ | $t_f = 1.7$ | $t_f = 1.8$ | $t_f = 1.9$ | $t_f = 2.0$ |  \n| -------- | ------- | ------- |  ------- |  ------- |  ------- |  ------- |  ------- |\n| 100  | $6\\times 10^{-4}$ | $1.04\\times 10^{-3}$ | $2.14\\times 10^{-3}$ |  $1.09\\times 10^{-2}$ | $3.41\\times 10^{-2}$ | $6.99\\times 10^{-2}$ | $1.09\\times 10^{-1}$ |  \n| 1000 | $8\\times 10^{-6}$ | $2.69\\times 10^{-4}$ | $1.74\\times 10^{-3}$ | $1.05\\times 10^{-2}$ | $3.32\\times 10^{-2}$ | $6.83\\times 10^{-2}$ | $1.06\\times 10^{-1}$ |  \n\nThese results shows that the extrapolation capacity of Neusa depends on the accuracy of the field in the training region. However, overfitting in the training region does not necessarily reduce the error in the more distant regions.\n\n### References\n\n[1] Brougois, A., Bourget, M., Lailly, P., Poulet, M., Ricarte, P., & Versteeg, R. (1990, May). **Marmousi, model and data**. In *EAEG workshop-practical aspects of seismic data inversion* (pp. cp-108). European Association of Geoscientists & Engineers.\n\n[2] Mishra, S., & Shekar, B. (2025). **Wave Equation and Neural Solvers: What Works and What Needs to Change?**. In 86th *EAGE Annual Conference & Exhibition* (Vol. 2025, No. 1, pp. 1-5). European Association of Geoscientists & Engineers.\n\n[3] Ding, Y., Chen, S., Miyake, H., & Li, X. (2025). **Physics-informed Neural Networks with Fourier Features for Seismic Wavefield Simulation in Time-Domain Nonsmooth Complex Media**. IEEE Transactions on Geoscience and Remote Sensing.\n\n[4] Zhongkai, Hao, et al. **\"Pinnacle: A comprehensive benchmark of physics-informed neural networks for solving PDEs.\"** Advances in Neural Information Processing Systems* 37 (2024): 76721-76774.\n\n[5] Liu, Dong C., and Jorge Nocedal. **\"On the limited memory BFGS method for large scale optimization.\"** *Mathematical programming* 45.1 (1989): 503-528."}, {"review_id": "0YnRrOCDoU", "reviewer": "Reviewer_4ytc", "summary": "The paper introduces NeuSA (Neuro-Spectral Architectures), a new class of Physics-Informed Neural Networks (PINNs) designed to solve complex partial differential equations (PDEs), including time-dependent, nonlinear, and non-homogeneous types. Traditional MLP-based PINNs often struggle with spectral bias (favoring low frequencies) and causality violations. NeuSA overcomes these issues by:\n\n- Projecting PDEs onto a spectral basis, capturing high-frequency dynamics.\n- Integrating with Neural ODEs, which preserve temporal causality.\n- Using classical spectral methods for initialization, starting training closer to the target solution.\n\nNeuSA demonstrates superior performance on various benchmark PDEs (e.g., wave equation, Burgers’ equation, sine-Gordon equation), achieving higher accuracy, faster convergence, and better temporal consistency than other PINNs.\n\nThe architecture naturally satisfies boundary and initial conditions without extra training cost. Future work may extend NeuSA to more complex PDEs like Navier-Stokes and explore alternative spectral bases and numerical solvers to further improve performance and generalization.", "questions": "The authors should justify their claims by testing NeuSa in more complex setups and on a wider variety of PDEs (high dimensions, complex transient dynamics, etc.)", "limitations": "The major limitation of this work is the lack of an extended benchmark that would convince the readership of the claimed advantages of NeuSA. These kinds of benchmarks are nowadays necessary in order to position a new PINN version within the wide PINN's community.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\nthe spectral formulation, alongside the introduction of NeuralODE, seem to increase the accuracy in standard benchmark and preserve causality.\n\nWeaknesses: \nNeuSA has not been tested in extreme case scenarios (very high dimensions, complex dynamics, etc.). It is therefore complex to assess NeuSA position within the vast PINN panorama.", "was_revised": true, "final_justification": "The authors' rebuttal clarified some aspects and showed further comparison with other existing models. However, most of those comparisons were limited to 1D and 2D cases. I encourage the authors to benchmark NeuSA against 3D challenging cases.", "rebuttal": "We thank the reviewer for their observations and feedback. It is encouraging for us that they appreciate the strength of our spectral formulation, and that they recognize the increase in accuracy in many problems, as well as the positive effects of our causal structure. \n\n## Questions\n\n> **The authors should justify their claims by testing NeuSa in more complex setups and on a wider variety of PDEs (high dimensions, complex transient dynamics, etc.)**\n\nAlthough we understand the reviewer’s concern in testing NeuSA in more complex scenarios, we believe the equations used in the paper already pose very significant challenges to current PINN architectures, as pointed out by Reviewer 3zAY. Likewise, we believe the weak performance of many of the reference architectures on the chosen benchmarks, which include reflections and shocks, demonstrates that the chosen benchmarks do in fact display complex dynamics.\n\nWe point out that these benchmarks go substantially beyond those used in the introduction papers for our baselines. These often deal with linear, translation-invariant problems, which NeuSA can solve at initialization to near-perfect accuracy by construction. Furthermore, most studies proposing new strategies for solving PINNs do not address high-dimensional time-dependent problems. This includes works such as QRes [1], sf-PINN [2], RoPINN [3], and PINNsFormer [4]. In particular, the PINNacle benchmark [5] notes that \"when systems have higher dimensions (3D Poisson equation) or longer time spans (2D Heat equation), all methods fail to solve, highlighting the difficulties associated with complex geometries and multi-scale systems.\" \n\nNevertheless, we believe the concern about the architecture's performance in high-dimensional problems is valid. To demonstrate NeuSA’s capabilities in higher dimensions, we have run an additional experiment on the wave equation in three spatial dimensions over a two-layered medium. The results can be found in the table below:\n\n| Model | Loss | rMAE | rMSE | Training Time (s) |\n|---|---|---|---|---|\n| PINN | $0.371$ | $0.145$ | $0.0613$ | $6168$ |\n| QRes | $0.078$ | $0.0298$ | $0.0167$ | $8758$ |\n| FLS | $0.313$ | $0.124$ | $0.0563$ | $6288$ |\n| NeuSA| $0.00107$ | $0.00561$ | $0.00435$ | $737$|\n\nAs we can see in the table above, Neuro-Spectral Architectures are able to achieve the smallest rMSE with a much smaller training time than any other reference architecture. This experiment and a longer discussion of its results will be included in the final version of the paper.\n\nWe also appreciate the reviewer’s comment regarding the need to assess NeuSA in a wider range of problems. We would like to emphasize, however, that the current benchmarks cover a relatively wide range of problems and in many ways improve upon basic test cases typically found in the PINNs literature. In particular, our numerical tests consist of:\n\n1. The 2D wave equation in heterogeneous media, which introduces discontinuities in wave speed and demands high-frequency resolution to capture reflections and mode conversions accurately. Our benchmarks address substantially more complex geometries when compared to those reported in [2,3,4], all of which focus on the 1D wave equation in a homogeneous medium. Our benchmark is also more complex than its analogue in PINNacle [5], due to the presence of reflections.\n\n2. The 1D sine-Gordon equation, which is characterized by strong nonlinearity and has proved to be a challenging test for the baseline architectures. This equation can be regarded as a more difficult, nonlinear counterpart to the 1D linear wave equation considered in many of our references.\n\n3. The 2D Burgers’ equation with periodic initial condition, which captures nonlinear advection and the evolution of complex flow structures. This generalizes the 1D problem as it appears in QRes [5], and is almost identical to the setup in PINNacle [4].\n\nThe examples provided in the paper were carefully selected to evaluate NeuSA's ability to overcome key limitations of standard PINNs, including spectral bias and poor performance in time-dependent problems. We have also focused on PDEs of practical interest with nonlinearities and inhomogeneous coefficients. Many of the other PDEs listed in the PINNacle benchmark [4] are linear and translation invariant; as mentioned above, NeuSA could solve them with zero training time, which we believe would be an unfair comparison.\n\nTo further demonstrate NeuSA’s performance on challenging benchmarks, we have run a second additional experiment on the 2-D wave equation, considering the Marmousi model [6]. This is a highly irregular acoustic model that induces highly complex reflections. This problem is a canonical benchmark in acoustic geophysics and is challenging for even specialized PINNs [7].\n\nWe tested NeuSA on the 2D wave equation using a cropped version of the Marmousi model. By doubling the\nfrequency content and training for 3000 epochs, we achieved an RMSE of 9.05e-2 with a training time of 923 seconds. This result is comparable to the framework tackled by [8], whose architecture is specialized to simulating seismic wavefields for nonsmooth complex media. This model has a much larger set of highly optimized hyperparameters and was trained using time-domain decomposition over 45000 epochs with 80000 collocation points, and requiring hours of training in our current setup. We therefore conclude that NeuSA has a strong performance in this scenario when compared to most of the literature.  We will include these results to the final version of our paper.\n\n### References\n\n[1] Bu, J., & Karpatne, A. (2021). **Quadratic residual networks: A new class of neural networks for solving forward and inverse problems in physics involving pdes.** In Proceedings of the 2021 SIAM International Conference on Data Mining (SDM) (pp. 675-683). Society for Industrial and Applied Mathematics.\n\n[2] Wong, J. C., Ooi, C. C., Gupta, A., & Ong, Y. -S. (2024) **Learning in Sinusoidal Spaces With Physics-Informed Neural Networks.** IEEE Transactions on Artificial Intelligence, vol. 5, no. 3, pp. 985-1000.\n\n[3] Wu, H., Luo, H., Ma, Y., Wang, J., & Long, M. (2024). **ROPINN: Region optimized physics-informed neural networks.** Advances in Neural Information Processing Systems, 37, 110494-110532.\n\n[4] Zhao, Z., Ding, X., & Prakash, B. A. (2023). **PINNsformer: A transformer-based framework for physics-informed neural networks.** arXiv preprint arXiv:2307.11833.\n\n[5] Zhongkai, H., Yao, J., Su, C., Su, H., Wang, Z., Lu, F., Xia, Z., Zhang, Y., Liu, S., Lu, L., & Zhu, J. (2024). **PINNacle: A comprehensive benchmark of physics-informed neural networks for solving pdes.** Advances in Neural Information Processing Systems, 37, 76721-76774.\n\n[6] Brougois, A., Bourget, M., Lailly, P., Poulet, M., Ricarte, P., & Versteeg, R. (1990). **Marmousi, model and data.** In EAEG workshop-practical aspects of seismic data inversion (pp. cp-108). European Association of Geoscientists & Engineers.\n\n[7] Mishra, S., & Shekar, B. (2025). **Wave Equation and Neural Solvers: What Works and What Needs to Change?.** In 86th EAGE Annual Conference & Exhibition (Vol. 2025, No. 1, pp. 1-5). European Association of Geoscientists & Engineers.\n\n[8] Ding, Y., Chen, S., Miyake, H., & Li, X. (2025). **Physics-informed Neural Networks with Fourier Features for Seismic Wavefield Simulation in Time-Domain Nonsmooth Complex Media.** IEEE Transactions on Geoscience and Remote Sensing."}, {"review_id": "qt5DWzbL6w", "reviewer": "Reviewer_aWbv", "summary": "This paper introduces a new method in the category of physics informed neural network (PINN) inspired by neural ode and spectral methods that can handle the causality problem faced by traditional PINNs. Instead of implementing the physical loss globally with time as an input variable, this method constructs a physical loss for each discretized step in time evolution, and expresses the result as an integral of time steps. The spatial part of the solution function is treated as a linear superposition of basis functions. In this case, Fourier basis are chosen as basis functions. Networks are applied to learn time-dependent coefficients at discretized time steps. This method shows better performance on three equations in comparison with other PINN approaches, and can generalize well to future steps unseen during training.", "questions": "1. Are the basis functions in use fixed to be Fourier basis, or trained as networks? In lines 132-133 of main text, the paper indicates that the basis is \"initialized\" as Fourier basis. In lines 151-152 of main text, the paper also mentions that the basis is \"trained\". However, in equations 9 and 10, the dependency on trained parameters $\\theta$ are not present in the basis functions $\\bf{b(x)}$, indicating that those functions are not trained. In equations 25 and 26 in appendix B.1.2, the derivatives of exponential functions are calculated, also indicating that the basis functions are fixed to be Fourier. Can authors give a consistent description on basis functions?\n2. The description on boundary condition fulfillment is confusing. In lines 132-133 of main text, the paper proposes to use basis functions to fulfill boundary conditions. However, boundary conditions are not satisfied with the choice of basis functions. For example, non-homogeneous Dirichlet boundary conditions should be satisfied by adding a specific term in the format of solution, so that the boundary conditions are homogenized. Can authors explain more on how boundary conditions are satisfied by the choice of basis functions?\n3. In equation 11 of main text, a wave equation is shown to be solved in $\\mathbb{R}^2 \\times [0, 2]$, which means in spatial dimensions, the domain expands to infinity. However, figure 8 of appendix, the spatial solution domain is shown to be $[-2, 2]^2$. If the spatial domain is extended to infinity, then the boundary condition is set to be 0 at infinity. If the spatial domain is now $[-2, 2]^2$, then what is the boundary condition imposed here?", "limitations": "The limitations of this work are discussed in the paper.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\t\n1. The proposed method shows good extrapolation ability to unseen steps in the future, indicating that the model learns real \"physics\".\n2. A good convergence speed is achieved by this method, in comparison with other methods mentioned in this paper.\n3. The causality problem is solved with an approach like Neural ODE.\n\t\nWeaknesses:\n\t\n1. Some parts of this paper on boundary conditions and basis functions are confusing and inconsistent. A clarification is required. Please see question section for the most concerning points for the reviewer. \n2. A Runge-Kutta method is required for integrating the steps and obtain the result, which might make the inference slow.", "was_revised": true, "final_justification": "The rebuttal from authors answer most of my questions. Therefore I raise 1 point in the rating to recommend acceptance.", "rebuttal": "We thank the reviewer for their thoughtful observations and detailed feedback on the paper. We appreciate their recognition of the method's solution to the causality problem and the method's ability to extrapolate to unseen times, as well as their highlight of our convergence speed.\n\n## Weaknesses\n\n> 1 - **Some parts of this paper on boundary conditions and basis functions are confusing and inconsistent. A clarification is required. Please see question section for the most concerning points for the reviewer.**\n\nWe address these concerns in more detail below, when responding to each question individually. In short, we recognize that the description of these elements could be improved and will deliver a more detailed explanation for the paper's final version.\n\n> 2 - **A Runge-Kutta method is required for integrating the steps and obtain the result, which might make the inference slow.**\n\nThis is a reasonable concern and connects to a question by reviewer 3zAY. We discuss this briefly in the paragraph starting at line 273:\n\n\"Although built on the inherently slower NODE framework, NeuSA trains significantly faster than MLP-based PINNs\"\n\nIn essence,  inference and training of these models is indeed slower, but the much faster convergence rate (in terms of number of epochs) generally more than compensates for this effect. For the paper's final version we will include an additional snippet expanding on this point, including the following table with the inference time (in seconds) for each architecture.\n\nWe compare the time needed for each architecture to predict the solution over a spacetime grid, using the same discretizations as in the experiments. The results are as follows:\n\n| Architecture | Wave 2D (101x101x201) | Burgers 2D (201x201x101) | Sine-Gordon (201x201 |\n| -------- | ------- | ------- | ------- |\n| NeuSA | 0.12 | 0.11| 0.10 |\n| MLP | 0.04 | 0.28 | 0.0005 |\n| QRes | 0.11 | 0.34 | 0.002 |\n| FLS | 0.04 | 0.28 | 0.0006 |\n| PINNsFormer | 1.29 | 2.57 | 0.02 |\n\nAs can be seen from the table, NeuSA's inference time is dominated by the numerical integration process, while the baseline architectures process every point in the spatio-temporal grid simultaneously. This results in slower inference when iterating over large amounts of points. Throughout our experiments, we mitigate this effect for the reference architectures by randomly sampling the domain in the experiments with the 2D Wave/2D Burgers equations.\n\n## Questions\n\n> 1 - **Are the basis functions in use fixed to be Fourier basis, or trained as networks? In lines 132-133 of main text, the paper indicates that the basis is \"initialized\" as Fourier basis. In lines 151-152 of main text, the paper also mentions that the basis is \"trained\". However, in equations 9 and 10, the dependency on trained parameters  are not present in the basis functions, indicating that those functions are not trained. In equations 25 and 26 in appendix B.1.2, the derivatives of exponential functions are calculated, also indicating that the basis functions are fixed to be Fourier. Can authors give a consistent description on basis functions?**\n\nWe agree with the reviewer that the discussion of classical bases and their trainable counterparts in the paper requires an improvement, and we will address this in our final revision.\n\nIn our experiments, we have focused on using the Fourier bases with fixed frequencies. Though at first we attempted to use trainable bases, a more detailed ablation showed this approach underperformed when compared with a fixed basis. We came to the conclusion that these tools demanded a more thorough theoretical and practical analysis, and will add this as a future direction of research.\n\n> 2 - **The description on boundary condition fulfillment is confusing. In lines 132-133 of main text, the paper proposes to use basis functions to fulfill boundary conditions. However, boundary conditions are not satisfied with the choice of basis functions. For example, non-homogeneous Dirichlet boundary conditions should be satisfied by adding a specific term in the format of solution, so that the boundary conditions are homogenized. Can authors explain more on how boundary conditions are satisfied by the choice of basis functions?**\n\nWe appreciate this remark. In our study, we restricted our attention to periodic and homogeneous Dirichlet and Neumann conditions, for which enforcing the boundary conditions through the choice of basis functions is straightforward. We will clarify that the boundary conditions are homogeneous in the final version. \n \nWe also note that when the domain geometry or boundary conditions cannot be enforced automatically by the choice of basis, NeuSA may still enforce these constraints by means of Physics-Informed loss terms, as used by normal PINN architectures. We also intend to emphasize this aspect more clearly.\n\n> 3 - **In equation 11 of main text, a wave equation is shown to be solved in $\\mathbb{R}^2 \\times [0,2]$, which means in spatial dimensions, the domain expands to infinity. However, figure 8 of appendix, the spatial solution domain is shown to be $[-2,2]^2$. If the spatial domain is extended to infinity, then the boundary condition is set to be 0 at infinity. If the spatial domain is now $[-2,2]^2$, then what is the boundary condition imposed here?**\n\nThis is a pertinent observation. The setup we are investigating here is very common in the seismics community, which is to approximate solutions to the wave equation over an infinite domain by solving it on a large domain with homogeneous Dirichlet or Neumann boundary conditions. This approach is valid as long as the domain is large enough so that the waves do not reach the boundary. Alternative approaches, including using absorbing boundary conditions, are an important direction of future research.\n\nWe will clarify this in the final version of the paper and explain more carefully the relationship between the experiment's domain and the original physical problem."}]

Author rebuttal:

Meta-review: This paper demonstrates several important strengths that justify its acceptance. Most notably, the proposed method exhibits strong extrapolation ability to unseen future steps, indicating that it successfully captures the underlying physical principles rather than merely fitting the training data. This is a critical feature for ensuring robustness and generalizability in scientific and engineering applications. In addition, the method achieves a faster convergence speed than other approaches discussed in the paper, highlighting its efficiency and practical utility, especially for large-scale problems where computational cost is significant. Furthermore, the authors effectively resolve the causality issue by employing a strategy inspired by Neural ODEs, thereby maintaining temporal consistency and physical plausibility in the learned dynamics. Equally important, the authors have carefully and appropriately addressed the reviewers’ questions, providing clarifications and revisions that further strengthen the manuscript. Taken together, these merits make the paper a valuable contribution, and I recommend acceptance.
Decision: Accept (poster)