Enhancing Physics-Informed Neural Networks Through Feature Engineering
Authors: Shaghayegh Fazliani
Zachary Frangella
Madeleine Udell
Venue: Accepted by TMLR ()

Abstract: Physics-Informed Neural Networks (PINNs) seek to solve partial differential equations (PDEs) with deep learning.  Mainstream approaches that deploy fully-connected multi-layer deep learning architectures require prolonged training to achieve moderate accuracy, while recent work on feature engineering allows higher accuracy and faster convergence. This paper introduces SAFE-NET, a Single-layered Adaptive Feature Engineering NETwork that improves errors with far fewer parameters than baseline feature engineering methods. SAFE-NET returns to basic ideas in machine learning, using Fourier features, 
a simplified single hidden layer network architecture, and an effective optimizer that improves the conditioning of the PINN optimization problem. Numerical results show that SAFE-NET converges faster and typically outperforms deeper networks and more complex architectures. It consistently uses fewer parameters --- on average, 53% fewer than the competing feature engineering methods and 70-100$\boldsymbol{\times}$ fewer than state-of-the-art large-scale architectures --- while achieving comparable accuracy in less than 30% of the training epochs. Moreover, each SAFE-NET epoch is 95% faster than those of competing feature engineering approaches. These findings challenge the prevailing belief that modern PINNs effectively learn relevant features and highlight the efficiency gains possible through feature engineering.

Published in Transactions on Machine Learning Research (11/2025)
Enhancing Physics-Informed Neural Networks Through Fea- ture Engineering
Shaghayegh Fazliani
fazliani@stanforddeee edu
Deportment of Mathematics, Stanford University, Stamford, CA, USA
Zachary Frangella
Strangella@sttnnorrdco
Department of Management Science & Engineering, Stanford Universiby, Stanford, CA, USA
Madeleine Udell
maddll@sttafordd.cm edu
Deportment of Management Science & Engineering, Stanford University, Stanford, CA, USA ICME, Stanford University, Stanford, CA, USA
Reviewed OpenReview: https: openrevvee.......... net/ formm??m?? id= Jd=J250
Abstract
Physics-Informe Neural Networks (PINNs) seek to solve partial differential equations (PDEs) with deep learning. Mainstream approaches that deploy fully-connected multi-layer deep learning architcctures require prolonged training achieve moderate accuracy, while recent work on feature engineering allows higher accuracy and faster convergence. This paper introduces SAFE-NES NEN,, Single-layered Adaptive Feature Engineering NETwork that improves errors with far fewer parameters than baseline feature engineering methods. SAFE- NET returns to basic ideas in machine learning, using Fourier features, a simplified single hidden layer network architecture, and an effective optimizer that improves the conditioning of the PINN optimization problem. Numerical results show that SAFE-NET converges faster and typically outperforms deeper networks and more complex architectures. It consistently fewer parameters on average, 53% fewer than the competing feature engineering methods and 70-100x fewer than state-of-the--at large-scale architectures while achieving comparable accuracy in less than 30% of the training epochs. Moreover, each SAFE-NET epoch is 95% faster than those of competing feature engineering approaches. These findings challenge the prevailing belief that modern PINNs effectively learn relevant features and highlight the efficiency gains possible through feature engineering.
1 Introduction
Partial Differential Equations (PDEs) underpin scientific modeling but remain notoriously challenging solve. Classical numerical methods struggle with high dimensionality and nonlinearity, while analytical solutions are rare. Traditional numerical methods are supported by software libraries such FEniCS (Alnaes et al. 2015). deal.1l (Arndt et al. 2021) and PETSe (Balay et al. 2025), but still considerable domain expertisc is required for effective implementation. More recently, LLM-bbsed approaches have been explored for aulomating PDE solver generation, combining symbolic reasoning and code synthesis to produce solvers Li et al. 2025 Soroco et al. 2025 Fazliani Udell 2025 Such approaches might struggle when the target PDE lacks prior numerical treatment or is underreppesentee in the training data, limiting their generality. Physics- Physics-Informed Neural Networks (PINNs) (Raissi et al. 2019a Karniadakis al.) [2021) are a promising alternative, leveraging neural networks to approximate PDE solutions through residual minimization. They aim to solve PDE systems of the form
Diu(x), = 0, B B[u(x), =0
x E Ω x E aQ
(1)

---

Published in Transactions on Machine Learning Research (11/2025)
where D represents the differential operator defining the PDE, B represents a general boundary condition operator, and Ω S IRd is the domain of the PDE. By avoiding mesh generation, PINNs offer flexibility for forward/inverse problems and high-dimensional settings.
Despite their potential, PINNs face a fundamental challenge- they are difficult to train Krishnapriyan et al. 2021 Rathore et al. 2024). Research shows that the differential operator the residual loss induces ill-conditioning (De Ryck et al. 2023; Rathore al.) 2024], leading to a poor optimization landscape and slow convergence for popular first-order optimizers such Adam (Kingma & Ba 2014]. Thus, successful PINN training can be time-consuming and fiddly, limiting the of PINNs.
Recent work has developed several strategies to improve PINN training, including feature engineering, Feature engineering endows the network with additional features that better caplure the inductive bias of the learning task. A range of feature engineering approaches, from Fourier features (RFF) (Wang et al. 2020 to radial basis function features (RBF) (Zeng al., 2024a), have been proposed. This focus on feature engincering represents a departure from the dominant trend in modern deep learning, where end-to-end learning has largely replaced hand-crafted features. However, the scientific computing domain presents unique challenges that may favor explicit feature design. Unlike image recognition natural language processing, where optimal feature representations are unknown a priori, PDE solutions have well-understood mathematical structures rooted in Fourier analysis, differential geometry, and physics principles. This domain knowledge provides strong theoretical guidance for feature design that is often absent in other ML applications. Recent work in domain-specific architectures from Graph Neural Networks for relational data Bronstein et al. (2017) to Transformer attention mechanisms for sequences Vaswani knowledge outperform generic architectures.
(0177) suggests that incorporating structural
However, prior work on feature engineering generally sulfers from one more of the following four limitations: 1) they impose rigid priors (e.g. the features are fixed or random functions), 2) they require hyperparameter tuning (e.g. determining kernel hyperparamcters in RBFs), 3) they are often computationally expensive, and 4) they fail to integrate domain knowledge such as boundary or initial conditions. Thus, while existing feature engineering techniques can improve performance under certain conditions, they can be PDF-specific, expensive, and sensitive to hyperparameters.
PINN
RBF-P
W-PINN
RBA-PINN
FLS-PINN
RBF
RFF
SAFE-NET
L2RE
10-1
10-2
Iterations
x10
Figure 1: L2RE for the PDE with SAFE-NET and basclines methods using Adam. For detailed nunerical results see Table 19 in Appendix
C.2
To address these shortcomings we introduce SAFE-NET, a feature engincering framework for PINNs that combines:
(1) Well-conditioned adaptive Fourier basis terms as features with trainable frequencies was A, and amplitudes adapt to PDE-specific dominant frequencies while mitigating spectral bias.
(2) Domain knowledge features derived from initial and boundary conditions (if available), providing added benefits in boosting accuracy.
2

---

Published in Transactions on Machine Learning Research (11/2025)
(3) Normalization of features prior to network imput stabilizes quasi-Newton optimizers (c.g., L-BFGS) and other advanced optimizers and prevents divergence. Many other methods introduce sensitive architectural or problem- hyperparameters or features, causing instability across some of our tested PDEs. In other words, SAFE-NET can safely be optimized by high-performance optimizers.
SAFE-NET does not diroctly target the optimization landscape like prior fcature engincering approaches such as RFFs or RBFs. Instead, it seeks to improve the inductive bias of the PINN by augmenting the initial data representation with well-conditioned Fourier features and domain knowledge. These design choices make SAFE-NET particularly effective for PDEs without shocks or discontinuities, where Fourier bases work best. Interestingly, this inductive bias also provides implicit preconditioning effect that leads to betlee-condiitonne optimization problem and facilitates training.
Contributions. We highlight the contributions of this paper:
Our new conputationally--ffcceent feature engineering method, SAFE-NET, offers better inductive bias than cxisting feature engincering methods.
We demonstrate empirically that SAFE-NET implicitly preconditions the loss landscape, leading to faster, more stable convergence.
SAFE-NET achieves runtime performance comparable to non-feature engincering methods with very low parameter count, while other feature engineering approaches are slower (Table [ and Figure 4
Experiments across wide variety of PDEs, including (i) linear, mildly nonlinear, and highly nonlinear dynamics, (ii) time-dependent and steady-state problems, and (fii) PDEs with predominantly smooth behavior well as discontinuous shock-forming PDFs, in dimensions < 3 and baseline methods (Table 2 show that on average, SAFE-NET yields the best or comparable performance.
Deep learning has been seen as a promising tool for solving PDEs. However, our work with SAFE-NET shows that traditional techniques like feature engineering can achieve lower error rates and faster training times than conventional multi-layer networks, particularly suited to PDEs whose solutions are sufficiently smooth. While recent. large-scale approaches such as PINNsFormer Zhao et al. (2024), PirateNets Wang et al. (2024), and operator learning approaches such FNO Li al. (2020b) and its variations may achieve superior accuracy in certain PDEs through complex architectures, these approaches
extremely expensive they use very large networks: e.g., 454k parameters for PINNsFormer Zhao et al. 2023). 720k parameters for PirateNets Wang et al. (2024). and approximately 1M parameters for FNO and similar approaches. This paper shows that for many PDE tasks, simple feature engineering can perform as well or better than complex architectures, challenging the preconception that deeper networks effectively learn important problem features.
2 Insights into Feature Engineering in PINNs
This section explores the two primary approaches to feature engineering, Fourier-based and non-Fourier feature mappings. To better understand the advantages SAFE-NET offers, we highlight the strengths, limitations, and applicability of these approaches different PDE classes.
2.1 Fourier-Based Feature Engineering
Fourier feature mappings leverage the spectral properties of PDE solutions to enhance high- high-frequency learning. They aim address spectral bias the tendency of neural networks to favor low-frequency functions by transforming input coordinates into
expressive reprosentation. The most prominent approach, RFF Wang el al. (2011) uses the feature mapping
7(0) = [cos(BB) sin)
(2)
with fixed Gaussian weights E Ra drawn from AN(0, o2). With both cosine and sine terms, this mapping projects inputs into high-dimensional space where periodic and high-frequency patterns are more easily
3

---

Published in Transactions on Machine Learning Research (11/2025)
captured. Theoretical insights from Tancik al. (2020b) show that Fourier features help neural networks learn high-frequency functions in low-dimensional domains, which are typical in PDE applications. However, RFF's reliance on fixed, random frequencies limits its adaptation to PDE-specific spectral properties. Moreover, the Gaussian initialization B may not align with the dominant frequencics the solution, leading to reducod performance.
Another baseline method we study, RBA-PINN from Auagnostopoulos al. (2023), also incorporates Fourier feature embeddings for enforcing periodic boundary conditions in some of its variants. While the main contribution of RBA-PINN is the residual-based attention weighting scheme, their ablation studies demonstrate that Fourier features provide the most significant performance improvement among all components tested. Specifically, for the Allen-Cahn equation, RBA alone achieves 3.16 x 10-3 relative L error, while RBArie improves this by two orders of magnitude. The authors explicitly acknowledge that the coupling of the RBA scheme along with the Fourier feature embedding the most important component achieving low L2RE. We therefore include RBA-PINN as feature engincering bascline (for periodic boundary conditions), recognizing that while RBA's 's innovation lies in adaptive weighting, its west-performing configuration on PDEs with periodic boundaries relies heavily on Fourier feature engineering for its effectiveness.
2.2 Non-Fourier Approaches
For PDEs with sharp gradients discontinuities, Fourier features may struggle due to the Gibbs phenomenon; see Zeng et al. (2024b). The Burgers PDE, with its sharp discontinuity at x = 0, is well-suited to observe this behavior. Comparing numerical results for the Burgers PDE across methods in Table 3 show that RFF performs comparably to methods without feature engineering but underperforms RBF by over an order of magnitude.
Zeng et al. (0024a) addresses the Gibbs phenomenon by using Radial Basis Functions (RBFs)
|| - c11
@RBF(I; exp
(3)
20
where 6 denotes the center and controls the kernel width. RBF expansions can better approximate local, abrupt changes but computationally intensive due to kernel regression requirements and less suited for periodic or high-frequency PDE solutions Tancik et al. (020bb). The RBF fcature mapping function is defined as:
1= wi (RBF(E;
4 2z)
(4)
E ÓRBF (z; C, o)
where Ys are the centers of the RBFs (trainable parameters) and {w;'s Y's are weights for the feature mapping layer.
Polynomial Enhancement. A key limitation of pure RBF approaches is their difficulty in representing global trends and polynomial components that frequently appear in PDE solutions. To address this, RBF-P (RBF with Polynomials) from Zeng et al. (2024b) augments the RBF feature set with polynomial terms, so the feature mapping function is adjusted
Σx w."h CRRF cil) + w(t
(xx)
ΣnT ÓRBF([X (aBF((x-ci) c:)
where P is the polynomial function. In the feature mapping layer, it can be represented as (RBp(r!)
fi
(RBF(rT")
1 : IN
x
Wm
:
W
JN
@RBP (rk)
ΦRBF
with r; = -- c.. The polynomial terms capture global struct ure while RBFs handle local variations. The weights serve as Lagrange multipliers, enabling constraints on the RBF coefficients in the parameter space. This hybrid approach combines the local approximation power of RBFs with the global trend-fitting capability of polynomials, often resulting in improved accuracy over pure RBF methods. We compare against both RBF and RBF-P in our study (see Appendix A.2 for details on the experimental set-ups for the baseline methods).
4

---

Published in Transactions on Machine Learning Research (11/2025)
2.3 Comparisons & Practical Considerations
The choice between Fourier and non-Fourier features depends PDE characteristics. Fourier features ideal for smooth, periodic solutions but struggle with discontinuities. RBFs handle sharp discontinuities better but are computationally costly. Experiments show that the polynomial terms in RBF-P are also more effective in nonlinear equations like the Burgers Equation and Navier-Stokes; see Table 3 for numerical comparisons. Despite the disadvantage of Fourier features these tasks, SAFE-NET offers a sensible compromise, allowing trainable frequency parameters and domain-inspired features to improve inductive bias with unified computationalll--ffeeeet design.
3 Methodology
We introduce SAFE-NET in this section. We begin with some motivation from Fourier analysis.
3.1 Theoretical Background
Let f(x) R.d - R be function defined on a d-dimensional domain. Under mild regularity condi- tions, f(x) can be reconstructed from its Fourier transform using the inverse Fourier transform f(x) S f(x)e2xxixx dn, where f(x) the Fourier transform of f at frequency &. To approximate S, we can focus the dominant frequencies with large If(r)). Summing over these dominant frequencies, we obtain (xx Ex dominant 2xixx We can express the Fourier transform f(x) through its real and imaginary components to rewrite the approximation as
f(x) ~
Σ
(A, cos .x) x) + B₁ sin (2Tx x)),
K dominant
where An and B^ are real-valued coefficients derived from fK)). Fourier basis elements are effective as features when they include the dominant frequencics x of f(x).
SAFE-NET
Feature Generator Module
Normalization
Lineer Layer
Tanh Activation
Linear Layer
(x, t)
Outputs
Feature Generator
8. ((::i((,,))o0(((.) & cos((A,))) sin(se)t) s,si((()) sin(w,t), d;cos(2,x cos(or;))
Domain Knowledge Festures (f ava availablej
Figure 2: Diagram showing how SAFE-NET works for 1D time-dependent PDE as example. The Feature Generator Module has trainable frequencies and coefficients for more elfective feat ure selection.
For PDE solution u(x,t)= 2 + R. the 2D Fourier transform and its inverse are given by
// u(x, t)e 2n drdtt
a(wa, Ax)
u(z,t) /// Úi(we de d 2ni(us dwydhee

---

Published in Transactions on Machine Learning Research (11/2025)
where (War Ae) are spatial and temporal frequencies. Expanding the complex exponential yields the tensor product basis
( a & 2
[cos(2mw)a)r) + isin (2Twaa)) @ (00s(2m)e)) + ssn (2m
For notational simplicity we absorb 2x into wz. d in what follows. This expansion produces four real-valued basis functions per frequency pair (War, Ae) as
0 a (x, 1) = cos(waar) cost (Aet)
(5)
Ae (x, t) sin(w (Aet)
(6)
0g A (x, f) = cos(w sin Act)
(7)
04 @ (x, () sin(w sin()et).
3.2 SAFE-NET
Motivated by the considerations of Section through learnable frequencies (we 1= and amplitudes (a (00,c(0 1 to estimate
3.1 SAFE-NET implements the parametric basis in equations
u t) [[a 0
(z, + bl 10 1o (x,
E=1
0 +
1 (x, + d
10 24"
(x,
(9)
where 0 {wx)) ,1 119) al4), b(0) cl0), d are trainable parameters. The explicit cross-frequency terms in equation 9 capture the tensor product structure of the 2D Fourier basis.
Domain Knowledge Features. The solution u(x, t) to a PDE often inherits structure from the domain geomelry, boundary conditions, and physical invariants. SAFE-NET explicitly encodes this domain knowledge through features (xxx() automatically derived from boundary conditions, initial conditions, and known solution patterns for each PDEs in the feature generator module. The systematic feature extraction follows these rules:
For Initial Conditions: Features extracted by analyzing the functional form of initial conditions u(x, 0) (z)). If uo(x) contains specific functional components (e.g., trigonometric, polynomial, or exponential terms), corresponding features are included. Therefore, if u(z, = Σ. axfx(xr) then (0 t) = (x)) are included as domain features.
For Boundary Conditions: If homogeneous Dirichlet conditions u(0, t) = u(L,t) 0 are present, features satisfying these conditions naturally (e.g., sin( (g=)) are prioritized by SAFE-NET's frequency initialization (Appendix A.1.4).
Domain Knowledge features for the PDEs considered in this paper appear in Appendix A.1 As shown in Figure 2 depending on the availability of domain information, these features are concatenated with the Fourier basis terms before normalization and linear projection. Ablation studies the role of these extra added features (see Appendix C.1.6) show that they result in improvements varying across different PDE tasks. Overall, they provide more accuracy boosting when the initial and boundary conditions carry more information of the final solution while being a lightweight supplementary component of SAFE-NET. SAFE-NET uses only one hidden layer: given an input (x, t) e x R, the SAFE-NET network computes
folæ 1) w a WWiOWWSAP..NEH (æ, ()) + b1) + b2,
(10)
with parameters 0 (WSAFE-NET: W1, bj w22 b2), nonlinearity tanh(-). and learnable feature mapping @f ((USAPENET,, (ar, f)). Complete architectural specifications, parameter initialization schemes, and training

---

Published in Transactions on Machine Learning Research (11/2025)
procedures for SAFE-NET are provided in Appendix A.1 To ensure robust experimental validation, we conducted comprehensive ablation studies examining key design choices including sensitivity to the number of Fourier features (Figure 9 and 18) activation function selection for the hidden layer (Figure 17 and Table 15
sensitivity to our initialization strategy for Fourier coefficients and frequencies (Tables 17 and 18). number of layers (Figure 19 and normalization (Figure 20 These ablation results, presented Appendix explaining the choice of each component.
By combining the generality of Fourier features with known solution characteristics, SAFE-NET accelerates convergence. Figure 3 compares SAFE-NET with the best competing PINN architecture for each PDE in Table 2
15
OguaCoometioo Eror SSAFE SAAFE-NET Erroi]
[
055
Peffrmance Lewel Competitive 10.5
Dehind (cl
Rec Cale
Bogger
PDES
Competitor Eroor Figure For each PDE, bars show 1og SAFE-NET Emor) comparing SAFF-NET to the best competitor baseline method (left) and median competitor baseline method (right). Colors represent performance levels, and rank annotations show SAFE-NET's position among all methods. SAFE-NET ranks first on 4/8 PDEs. It achieves top-2 performance on 5/8 PDEs; the PDEs which SAFE-NET falls behind (3rd or 4th place) are Allen-Cahn, Burgers and Navier-Stokes, which are non-lincar and not shock-free.
Cost. Table and Figure illustaate the parameter count and runtime (per epoch) for different bascline methods. The setups used for each method in our experiments appear in Appendix C suffices for SAFE-NET, reducing its parameter count and improving speed.
One hidden layer
Table 1: Parameter count comparison between baseline methods and SAFE-NET. Using the same number of features (128), SAFF-NET achieves a significantly lower parameter count than competing feature engineering methods while remaining comparable to non-feature engineering approaches.
PINN
FLS-PINN
W-PINN RBA-PINN
5.3k
5.3k
5.3k
16k
RFF
RBF
RBF-P
SAFE-NET
14.5k
14.2k
14.7k
6.8k
Scalability. The number of additional features in SAFE-NET increases exponentially with the spatial dimension of the PDE. Hence SAFE-NET offers improves accuracy and computational advantages for lower- dimensional problems, where better coverage of frequency space improves approximation quality, but is not suited problems with high dimensions. Developing effective feature engineering methods for higher dimensional PINNS is important challenge for future work.

---

Published in Transactions on Machine Learning Research (11/2025)
10a
PANN
RBF-P
WAPNN
RBA
10a
Time (s)
10-1
10-2
0.2
0.4
0.6
0.8
1.0
Number of Training Samples
x10
Figure 4: Average runtime comparison across varying numbers of training samples for bascline methods and SAFF-NET using identical computational resources and no concurrent processes. SAFF-NET demonstrates efficiency over competing feature engineering methods, maintaining runtime close to the less computationally expensive non-feature engineering methods.
4
Related work
Much work has been done to improve the training of PINNs that take different approaches from feature engineering. Broadly, these approaches can be divided into three categories: architectural modifications, loss-reweighting, and optimizer design. Recent efforts to improve PINN accuracy and speed include:
Architectural Modifications. One way to improve the training procedure for PINNs is to modify the network architecture, improving the optimization landscape relative to the basic PINN, making training easier. Examples of this approach include the adapative activation functions of Jagtap et al. 2020 (A-PINNs), and specialized architectures designed to mitigale spectral bias (Li et al. 2020a)..
Loss Re-weighting. Another popular technique is loss reweighting. For certain PDEs, the residual loss tends to dominate the boundary loss in that the optimizer focuses too much on minimizing the residual loss, leading to solution that fails to satisfy the boundary conditions. To address this, techniques like W-PINN Wang et al. (2022) balance loss components through heuristic or learned weights to down-weight the residual loss and better fit the boundary loss.
Optimizer Design. Another popular approach to PINN training is develop more sophisticated optimizers that are more robust to ill-conditioning. Several notable proposals in this area are the natural gradients method of Müller & Zeinhofer (202), MultiAdam Yao et al. (2023), and NysNewton-CGG Rathore et al. 2024 These methods larget. ill-conditioning directly and use curvature information from the loss or the model to precondition the gradient. This leads to an improved optimization landscape locally, enabling the optimizer take better steps and progress faster. In contrast, SAFE-NET enjoys an implicit proconditioning effect that is global by incorporating trainable feature layer, SAFE-NET changes the PINN objective, globally changing the optimization landscape. The results in Section 5.3 show SAFE-NET enjoys a significantly better- conditioned optimization landscape. Thus, SAFE-NET can be further combined with more sophisticated optimization schemes to obtain further improvements.
5 Results
We provide overview of the experimental setup and present results across multiple benchmarks.
8

---

Published in Transactions on Machine Learning Research (11/2025)
5.1 Baselines
We compare two groups of PINN baselines: (i) without feature engineering (PINN, FLS-PINN, W-PINN) and (ii) with feature enginoering (RFF, RBF, RBF-P, RBA-PINN). Details on each bascline could be found in Appendix A.2
As noted earlier in Section large models such as PINNsFormer (454k parameters), PirateNets (720k+ parameters), operator learning methods (1M+ parameters), and several other PINN variants with hundreds of thousands of parameters demonstrate the performance gains possible with large and complex architectures. In contrast, our experiments focus on well-performing baseline methods with simple structures and low parameter counts (at or below 16k; see Table [ especially emphasizing feature engineering approaches that offer computational efficiency. This scope allows us to demonstrate that simple feature engineering also yields accurate solutions at a fraction of the computational cost of the largest models.
We experiment on eight PDE benchmarks in this section provided in Table 2 Depending on availability, the datasets for these tasks are either from PDE benchmarks such PDEBench Takamoto et al. (2024) and
PINNacle Hao et al. (2023) or implemented directly if unavailable online. More details on each PDE and its
source are provided in Appendix B An additional experiment on a synthesized non-homogeneous boundary
heat. task is included in Appendix C.3
Table 2: Overview of the tested PDEs. Following the conventions of PINNs literature, the dimensions provided in this table are spatial dimensions. An extra dimension should be added for time-dependent PDEs reach the correct number of variables.
PDE
Dimensions
Type
State
Wave
1D
Linear
Time-dependent
Reaction
1D
Nonlinear
Time-dependent
Diffusion
1D
Linear
Time-dependent
Heat
Linear 2D
Time-dependent
Convection
Linear
-eeendedn
Burgers
1D
Nonlinear
T'iee-deeende
Allen-Cahn
Nonlinear
Time-dependent
Navier-Stokes
2D
Nonlinear
Steady-state
5.2 Experiments
Following prior work (such as Rathore et al. (2024) and Zeng et al. (2024b)), ) we use a two-stage Adam L-BFGS optimizer for the first set of experiments. Recent loss- studies on PINNs show that first-order methods such as Adam make rapid progress away from saddle points but eventually slow down in the illcoonditioned valleys induced by differential operators, whereas quasi-Newton methods like L-BFGS improve local conditioning yet are attracted to saddles and often stall when started from scratch. Combining them inherits the best of both worlds: an Adam "warm-up" steers the parameters clear of saddles; switching I.BFGS then preconditions the Hessian and yields much faster local convergence. Nevertheless, pure L-BFGS can terminate carly bocause its strong- Wolfe line scarch sometimes returns a step size of zero, especially on stiff PDEs Our schedule therefore runs Adam first until progress slows down significantly and switches to L-BEGS until the maximum number of iterations is reached. A more detailed theoretical justification of the Adam + L-BFGS schedule appears in Rathore et al. (004))
Optimization Schedule (1). The first set of experiments uses Adam for either 10k, 20k, or 30k iterations, followed by L-BEGS for the remainder of the training run (total 40k). We found that for our tested set of experiments, 30k iterations of Adam followed by 10k itcrations of L-BFGS performs best. Adam has an initial learning rate of 0.001, decaying exponentially by 0.9 every 2k iterations, followed by full-batch L-BFGS until convergence (tolerance 10-10) if the maximum number of iterations is reached.
9

---

Published in Transactions on Machine Learning Research (11/2025)
Optimizer settings are identical across all PDEs and architectures; each experiment is repeated with five random initializations, and the model performance is assessed using the relative & error (I.2RE), defined as
1y Ull |yyl
V (1i 0:)
VE=1
where & is the predicted solution and is the true solution. Notably, while we run every method for the same number of epochs, SAFE-NET achieves faster per-epoch computation times than compeling approaches, reducing its total training time.
Table 3: PDE benchmark results comparing several PINN architectures in relative T.2 error. The best results are shown in Bold. However, numerically close results (i.e., within 30% of the best result, including the best result in Bold itself) are also highlighted in Blue This accounts for potential variance in implementation and hyperparameter sensitivity. Cells with (*) indicate L-BFGS divergence in our experiments.
PDE
PINN
FLS-PINN
W-PINN
RFF
Wave
6.62e-2
4.37e-2
3
7.39e-3
Reaction
4.98e-2 2
1.27e-1
4.52c-2
1.20e-2
Diffusion
9.65e 3
5.89e-2 2
1.43e
2.56e
Heat
6.34e-3
6.98e-7
9.98e
8.56e-3 3
Convection
1.96e 2
4.12c-2 2
9.77e
Allen-Cahn
4.98e-1
1.46e-0
1.29e-1
1.03e-2
Burgers
1.01e-2 2
6.16e-2
Navier-Stokes
7.49e-1
7.91e-1
6.76e-1 1
PDE
RBA-PINN
RBF
RBF-P
SAFE-NET
Wave
1.57-3
2.48e-2
2.13e-2
1.21e
Reaction
3.31e-2
1.98e 2
1.37e-2
9.93e
Diffusion
1.29e-4
1.02e
9.91e-
1.21e-1
Heat
7.35e-4
3.65e-3
7.21e-4
5.31e-4
Convection
8.32e-3
7.72e-2
7.02e-3
-.37e-3
Allen-Cahn
8.21u-8
1.06e
1.97e-3
Burgers
4.17e-4
3.95e-3
3.47c-4 4
2.67e-3
Navier-Stokes
4.86e-1
2.56e-1
2.98e-1
5.26e-1
Even after warming with Adam, L-BFGS can still cause divergence or instability on some PDEs and baseline methods. Corresponding entries are denoted by *** in the result table. There could be several reasons why this is possible. For instance, RFF uses fixed random frequencics drawn from a Gaussian distribution, which can lead very poor conditioning depending on the random initialization. The mismatth between random frequencies and the true solution frequencies can create ill-conditioned optimization landscapes that cause L-BFGS to diverge. Alternatively, in the case of W-PINN and Navier-Stokes, W-PINN's adaptive weighting scheme computes weights based NTK matrix traces (see Wang et al. (022) for the details or the computation). Due to the complex flow physics of Navier-Stokes, the NTK matrices could become ill-conditioned near boundary layers or rapid changes in the adaptive weights can destabilize L-BFGS.
SAFE-NET, on the other hand, is compatible with quasi-Newton methods such as L-BFGS: it can safely be optimized by them. Table 3 demonstrates that for shock-free PDEs (i.e., not Burgers or Navier-Stokes, and Allen-Cahn to some extent), SAFE-NET outperforms competing methods in L2RE by in 4/5 tasks, following closely behind RBF-P for the diffusion task. (See Figure g for visual demonstration). Allen-Cahn is *less nonlinear' than Burgers or Navier-Stokes, but also develop steep gradients and interfaces, especially in phase separation dynamics. In practice, these sharp transition layers are not true shocks in the hyperbolic PDE sense, but numerically they still present challenges similar to shocks including Gibbs phenomena for Fourier features in particular. Despite that, SAFE-NET manages to come third best after RBA-PINN and RBF-P for both Allen-Cate Cahn and Burgers tasks. This outcome is not surprising as SAFE-NET is not designed for PDEs with shocks. Particularly, Zeng al. (2024b) suggests that the RBF kernel handles the discontinuity at 0 in the Burgers equation more effectively.
10

---

Published in Transactions on Machine Learning Research (11/2025)
Experiments with Alternative Combinations. While the Adam + L-BFGS combination in Optimization Schedule (1) demonstrates superior performance compared to individual optimizers (see Appendix AA3)) our analysis reveals that L-BFGS consistently stalls before reaching the maximum iteration limit. To further improve SAFE-NET's performance, we developed optimizer schedules that restart L-BFGS after it stalls. The Tollowing combination, which we call (Adam + T-BFGS)², particularly successful choice for optimizing the SAFE-NET architecture. Hyperparameter analysis Appendix A.3 validates our optimizer choice.
Optimization Schedule (2). The experiments with (Adam L-BFGS)² also rum for 40k iterations. First, Adam is used warm--u up phase (Adam(1)), followed by full-batch L-BFGS until it stalls (which is always the case in our experiments) or hits the maximum number of iterations. Then Adam (Adam(2)) runs again followed by a final round of L-BFGS. Testing multiple switch points for both Adam(1) and Adam((2 suggest. that for (Adam + L-BEGS) rumning a shorter first warm-up phase 3k iterations of Adam, followed by L-BFGS until stalls (which happen within 10k iterations), then reverting to Adam and switching L-BFGS for the last epochs works best for the tested PDEs. Since Adam (1) runs a short warm-up phase, we grid-searched the learning rate y) E {10-1, [10-2,10-2, 10-3, 10-3, 10-1, 10-1,10-5) per PDE. However, for Adam(2), the learning rate schedule matches choice in Optimization Schedule (1), which starts with initial learning rate of 2.2 decaying exponentially by 0.9 every 2k iterations. Both rounds of L-BFGS in (Adam + L- L-BFGS)² match our choices for L-BFGS in Optimization Schedule (1).
10
Adam
L-BFGS
10
10-1
Gradient Norm
10-4
10-2
10-5
10-3
10-6
10
10
Iteration
x10
Ieeation
x104
(a) Diffusion
100
10
Adam L-BFGS
10-2
10-1
Gradient Norm
10-3
10-2
10
10-3
10-6
2
Iteration
x10*
Iteration
x10
(b) Allen-Cahh
Figure 5: L2RE and gradient norm for the diffusion and Allen-Cahn PDEs using SAFE-NET with (Adam + L-BFGS)2.
Figure demonstrates the performance of SAFE-NET + + the diffusion and Allen-Cahi PDEs, presenting the gradient norm plots as well. Table summarizes these numerical results, showing significant improvements on more challenging PDEs containing steep interface transitions, shocks, and discontimuities (Allen-Cahn) Cahn, Burgers, and Navier-Stokes) along with improvements on other PDEs matching the top results of Table 3 These results demonstrate the effectiveness of (Adam + L-BFGS)² when paired with SAFE-NET, leading to significant improvements compared with Adam L-BFGS, especially in the diffusion and Allen-Cahn tasks. With this improved optimizer, SAFE-NET beats all compared architectures
11

---

Published in Transactions on Machine Learning Research (11/2025)
TngurCCopettionn Erro SSFE NET Errori
[1]
LG-
Perormance t-3x]
Competiive (055-20)
Dehlidd (40
PDF
Figure 6: Comparison of SAFE-NET (Adam L-BFGS)² results for each PDE against the best baseline performance achieved using either (Adam L-BFGS) or (Adam + L-BFGS)2. This comparison selects the optimal result for each competing method across both optimization strategies, highlights SAFF-NET's consistent advantages on shock-free PDEs (now ranking first on 5/8 tasks) and validates the effectiveness of (Adam + I-BFGS)² with SAFE-NET. Compare with Figure
on 5/8 of the tested PDEs as well as showing improvements for all the nonlinear PDEs with shocks; see Figure for visual demonstration.
Table 4: L2RE for SAFE-NET with Optimization Schedule (1) and Optimization Schedule (2)
PDE
Adam + L-BFGS
(Adam + L-BFGS)²
Wave
1.21e-3
9.81e 4
Reaction
9
8.91e-3
Diffusion
1.21e-4
8.23e-5
Heat
5.31e-4
3.60e-4 -
Convection
4.37c-3 3
3.84c 3
Allen-Calm
9.97e-4
4.13e
Burgers
-.
9.87e-4
NS
5.26c-1
3.7le-1
Tables 5 and 6 appearing in Appendix A.3 contain complete numerical results for Optimization Schedule (1) and (2) along with other combinations tested. Notably, (Adam + T-BFGS)) fails to improve (and often degrades) performance for the majority of the baseline methods across different PDEs, making it unsuitable to compare PINN architectures. Figure S shows the average percentage of improvement or deterioration (in 12RF) for each method after switching from Optimization Schedule (1) to Optimization Schedule (2). To see an example of the loss curves, Figure compares the performances of the top 3 best performing methods on the diffusion PDE using Adam L-BFGS (Optimization Schedule (1)) and (Adam L-BFGS)² (Optimization Schedule (2)).
Finally, Figure g shows the effects of adding Fourier fcatures on the wave, heat, and Navicr-Stokes PDEs; the error gradually decreases as the number of SAFE-NET features increases until it reaches saturation point (at around 120-140 features for these problems). This observation (and similar observations in Appendix C.1 suggests that there no advantage adding additional fcatures after a certain complexity is reachod: they do not improve performance, but only increase the number of trainable parameters and computational cost. Detailed numerical results for other PDEs provided in Table 16
12

---

Published in Transactions on Machine Learning Research (11/2025)
Adam L-BFGS
(Adam L-BFGSY
109
10a
REA-PINN
RBF-P SATE-NET
10
10-1
L2RE
10
LZRE
10-3
10-2
10
10
Iterations
x10
lterations
Figure 7: L2RE for RBA-PINN, RBF-P, and SAFE-NET on the diffusion PDE using Adam L-BFGS (left) and AAdam+La L-BFGS)² (right). SAFE-NET, previously ranking second in performance on the diffusion task, is now outperforming the other methods using (Adam + L-BFGS)
Accrage Porcentoge inorormmennt (5)
22
25
167%
3.4%
4..5
Romel
Method
Figure 8: Average percentage of improvement deterioration (in L2RE) for each method after switching from Optimization Schedule (1) to Optimization Schedule (2). Figure 12 in the appendix shows percentages for every method on every PDE, demonstrating that SAFE-NET is the only method that produces results in consistent improvement on every PDE after switching to (Adam + IL-BFGS)2.
10a
10-1
Heat Diffusion
12RE
10.2
Navier-Stokes
10-3
10
50
100
150
200
250
Number of Features
Figure 9: L2RE decreases as number of features increases.
13

---

Published in Transactions on Machine Learning Research (11/2025)
5.3 Feature Engineering & Spectral Density Analysis
SAFE-NET significantly improves the conditioning of all tested PDEs. We analyze conditioning both at initialization (3k epochs) and after training for while (100k iterations) using SAFE-NET and PINN with Adam as the optimizer with an initial learning rate of 0.001 and an exponential decay rate of 0.9 every 2k steps (Experimental setup of Figure
This experiment uses Adam only since L-BEGS has preconditioning effects that would confound our attempt to understand the effect of architecture on problem conditioning. Figure 10 illustrates the spectral density on the Hessian for the wave PDE as an example, plotted from the same experiment as in Figure [ Analogous plots for other PDEs provided in Appendix C alongside spectral density calculation details.
We see that SAFE-NET shifts most cigenvalues values close to at initialization (Figure 10 top). Post-training spectral density plots (Figure 10 bottom) reveal dramatic conditioning improvements: the top eigenvalues for the wave problem decrease by a factor of 10*, while those for the heat and Burgers problems decrease by 102 (Figures 23 and 26 Further experiments (Appendix () show that SAFE-NET reduces both the number and density of large cigenvalucs across all problems and all bascline methods (scc Table 21
102
SAFE-NET
10-1
Density
10
10-7
10-10
10
10
102
103
104
Eigenvalue
(a) Early into Training
102
10-1
Density
10-4
10-7
10-10
100
10
102
103
104
10
106
Dignvvauue
(b) End of Training
Figure 10: Spectral density for the wave PDE using SAFE-NET and PINN at the early stages of training and at the end of training.
This improvement. in conditioning is a key factor explaining why SAFE-NET features are easier to optimize. Figure II compares SAFE-NET and PINNs performance on different loss components for the diffusion and Allen-Cahn PDEs (from the same experiment as Figure 5) as an example. Rathore et al. (2024) argues the residual loss is the main cause of the ill-conditioned loss landscape of PINNs. SAFE-NET significantly improves the conditioning of each loss component, including the residual loss.
14

---

Published in Transactions on Machine Learning Research (11/2025)
Loss Components (SAFE-NET)
Loss Components (PINN)
10
10
Li
10
Lbc
100
10
Lees
10-1
10-2
10-2-
10-3
10-3
10
10-4
10-5
10-5
Eeration
x10
Reration
X10*
(a)
Diffusion
(SAFE-NET)
Loss (PINN)
10
10
100
109
Lbc
10-1
Lres
10-1
Error
10
10-2
10
10-3
10
10-4
10-5
10-5
Reeation
x10
keration
x10
(b) Allen-Cahn
Figure 11: Comparison of loss components of Figure 5 with SAFE-NET and PINN with (Adam + L-BFGS)2.
SAFE-NET significantly improves each loss component for both PDEs.
6 Implicit preconditioning in SAFE-NET
Empirical results in Section 5.3 show SAFE-NET improves the conditioning of each problem both at the beginning and end of training. To develop better intuition into this empirical observation for the early stages of training, we consider a simple didactic sctting similar to Wang et al. 20200 where the network given by
u(x, 1) = w T ((t,1),
(11)
here (2. E R and ø is the fcature map.
We begin with a definition:
Definition 6.1. For neural network fo(z) with parameters 0 e RP, the tangent kernel θy : x Ω + R given by
(z' Vof(z') Vof(z))VVoo() fz(z))
Given an input dataset X E Radd the tangent kernel matrix is the matrix with entries
(0y(0)), Vel(x;)
i.j = 1,..,,..,1
where and zy are ith and jth rows of X.
The neural tangent kernel 07 e guxn is the fixed kernel defined as
( x) lim ((O)(2', x)], p-000
where p the number of trainable parameters and the expectation taken over the weights at initialization Jacot et al. (2018); Liu al. (2200). The neural tangent kernel matrix θ is defined analogously to the Langent kernel matrix.
15

---

Published in Transactions on Machine Learning Research (11/2025)
In the limit p - oo, neural net training with fo is equivalent to kernel regression with the NTK matrix OF As training is reduced to kernel regression problem, the convergence speed of gradient-based optimizers is controlled by the conditioning of the NTK matrix @ Jacot et al. (0118); ;; Liu et al. (2022). Thus, a better-conditioned NTK yields a better optimization landscape and faster convergence.
In the context of 11 and under idealized conditions, we continue the following calculations as an example build intuition into the empirical results suggesting that SAFE-NET features lead to a better conditioned NTK and, so, better optimization landscape. When and p are large, the spectrum of 0 is closely related the spectrum of the integral operator Te, (g) , obtain control over the Te, spectrum for given by
<))(dr) Wang et al. 2020). Thus, we shall
We begin by writing the tangent kernel function corresponding to 11 As 0 w, Vou(z, = o(z,t), therefore the tangent kernel function is given by
On ((x, 1) (a', ((x,t) o( ,t).
Let us focus on the Fourier basis features. Using the notation of Section 3.2
we have
N
(x =
10
(x,
j=1
where 1Ø! fj-1 are defined in equations 5 and l = I,....N N are the number of cach type of product
features. y
The kernel is sum of separable lunctions as
(,((
ΣΣ
(2', (t',t).
(12)
The eigenvalues 3 and eigenfunctions (xx t) satisfy
/
(z Maba (1) de' = Brb(x,)) 1)-
Substituting (n((2, 1), (x', from equation 12 we get
N
(2,t) / (2t,)))))) ())(x,t)
(13)
Set (0 ()/d)) dd dt' by the domain compactness assumption, so for + 0, we get
1 (x, 1)
ΣΣ ( (3 (z,
(14)
3
To calculate the eigenvalues, substitute equation
14 into equation
13 to get
(22)2
i. j.m
3
Now, 10 (5) forms an orthonormal basis, as each @ is chosen from an orthogonal Fourier basis and has unit norm. Thus,
4 ΣΣ
10 (3))
ΣΣ 19 66
a
16

---

Thus, the Published in Transactions on Machine Learning Research (11/2025)
From this display, we conclude that for each pair (i,)), af (69))= eigenvalues of Te. are 0 or equal (2) with corresponding eigenfunction (t () = aj o5 (x,4). As the non-zero eigenvalues correspond to the directions relevant to learning, we focus on them. Recall the (c))) the trainable amplitudes and are set to at initialization. Hence, at initialization and early stages, we expect the condition number of the matrix i. to be approximately 1. If u were in the infinite width limit, this would imply @% well-conditioned and fast convergence of gradient- based optimizers. Thus, this idcalized example shows the features selected by SAFE-NET can lead to better conditioning and faster convergence. In practice, finite network width and perturbations from domain knowledge features and normalization are expected to introduce some spreading in the eigenvalue distribution. However, experiments show that SAFE-NET maintains denser eigenvalue distribution around the theoretical prediction from the example (see Figure 10fa)), preserving the core conditioning benefits anticipated by the idcalized example and its insights. Early on in training, Figure 10 shows the eigenvalue distribution is relatively uniform. During training, SAFE-NET gradually adjusts the c2 c balancing different frequencies rather than letting any single Fourier mode dominate. As result, the cigenvalue distribution shifts outwards relatively slowly. Thus, even at the end of the training, Figure IO shows the eigenvalue distribution has not changed much from its initialization, so the landscape remains well-conditiondd.
3 3= (c 2
(2)
} are
7 Discussion and Future Work
Our results encourage a fresh look at the importance of feature engineering for PDEs. While machine learning trends have favored complex architectures, our work suggests that enginecred features can offer a useful implicit bias that deep architectures struggle to replicate. PDEs present unique challenges, including ill-conditioned solution spaces, where deep learning techniques tend to struggle. While newer architectures could offer improvements, our thcory shows engineered features can always improve problem conditioning. Hence we expect feature engineering will have lasting importance for solving PDEs, complementing advances in deep learning.
Fourier features align best with low-to-mooderute--dineiiii PDEs demonstrating smoother behavior; for strongly nonlinear or Mhock-forming problems, non-Fourier features (e.g., RBF-P) are often more suitable, consistent with our results. While SAFE-NET demonstrates significant improvements conditioning and convergence for PDEs whose solutions exhibit predominantly smooth behavior (including piecewise-smooth profiles without strong shocks), several avenues remain to improve feature engineering in PINNs. The first is multi-stage feature learning, following Wang & Lai (002)), multi-stage approach that further boosts performance by itcratively learning best fit on the residual error, using SAFE-NET as a basc learner could be developed. Another direction is better incorporating physical priors. Enforcing physical laws such as conservation principles and symmetries into the feature engineering process could improve the performance of SAFE-NET on complex problems. Features based on concepts like Noether's theorem or Hamiltonian mechanics could provide stronger inductive biases. Finally, would be interesting to also consider non-Fourier features like radial basis functions for handling sharp gradients, which could allow SAFE-NET to perform well on wider range of problems. The challenge here lies in maintaining numerical stability when combining these different types of features. Most importantly, extending feature engineering to truly high-dimensional problems is important future work, but it presents very different challenges that are currently beyond the reach of PINN-style mothods.
17

---

Published in Transactions on Machine Learning Research (11/2025)
References
M. S. Alnacs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The FEniCS Project Version 1.5. Archive of Numerical Software, 3, 2015. doi: 10.11588/anss aan.2005.100.20555.
Sokratis J. Anagnostopoulos, Juan Diego Toscano, Nikolaos Stergiopulos, and George Em Karniadakis. Residual-based attention and connection to information bottleneck theory in pinns, 2023, URL https:
//axxiv.org/as//2200 003399
Daniel Arndt, Wolfgang Bangerth, Denis Davydov, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Jean-Paul Pelteret, Bruno Turcksin, and David Wells. The deal.ii finite element library: Design, features, and insights. Computers & Mathematics with Applications, 81:407 422, January 2021. ISSN 0898-1221. doi:
200.m....02000.0222. 02..0.... URL http:/tp:///x..o doi 100g///... 1016/j camwa. 2020. .2..0.2.000 022
Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil M. Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Picrre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman Mclnnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose F. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang, PETSc Web page. hntps::/ org/ 2025. URL https://petscc.re org/
Michael M. Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geometric deep learning: Going beyond euclidean data. IEEE Signal Processing Magazine, 34(4):18 42, July 2017. ISSN 1558-0792. doi: 10. (0.1100/msp.2017.26889
URL http://dx doi org/10. 10. 1109/MSP 2017. 2693418
Tim De Ryck, Florent Bonnet, Siddhartha Mishra, and Emmanuel de Bézenac. An operator preconditioning perspective training in informe- machine learning. arXiv preprint arXtv:2310.0580) 2023.
Shaghayegh Fazliani and Madeleine Udell. Pde-sharp: Pde solver hybrids through analysis and refinement passes, 2025. URL
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 2010.
Gene H Golub and Gérard Meurant. Matrices, moments and quadrature with applications, volume 30. Princeton University Press, 2009.
Zhongkai Hao, Jiachen Yao, Chang Su, Hang Su, Ziao Wang, Fanzhi Lu, Zeyu Xia, Yichi Zhang, Songming Liu, Lu Lu, and Jun Zhu. PINNacle: A Comprehensive Benchmark of Physics-Informa Informed Neural Networks for Solving PDEs. aixiv preprint arXiv:2306. 08827, 2023.
Mare Härkönen, Markus
and Bogdan Raitä. Gaussian process priors for systems of lincar partial differential equations with constant coefficients, 2023. URL
Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. Advances in neural information processing systems, 31, 2018.
Ameya D. Jagtap, Kenji Kawaguchi, and George Em Karniadakis. Adaptive activation functions accelerate convergence in deep and plysics-informed neural networks. Journal of Compulational Physics, 404:109136, 2020.
George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Lin Yang. Physics- informed machine learning, Nature Reviews Physics, 3(6):422 440, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. Xi preprint 1412:111.2.00 6980, 2014.
18

---

Published in Transactions on Machine Learning Research (11/2025)
Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. In Advances in Neural Information Processing Systems, 2021.
Ke Li, Kejun Tang, Tianfan Wu, and Qifeng Liao. D3M: A Deep Domain Decomposition Method for Partial Differential Equations. IEEE Access, 8:5283-555 5294, 2020a.
Shanda Li, Tanya Marwah, Junhong Shen, Weiwei Sun, Andrej Risteski, Yiming Yang, and Ameet Talwalkar. Codepde: An inference framework for Ilm-driven pdc solver generation, 2025. URL https://arxiv org/ abs/2505.005 .07883
Zongyi Li, Nikola Kovachki, Kamyar Aziwadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint Xiv:20000. 08899, 2020b.
Chaoyue Liu, Libin Zhu, and Misha Belkin. On the linearity of large non-linear models: when and why the tangent kernel is constant. Advances in Neural Information Processing Systems, 2020.
Chaoyue Liu, Libin Zhu, and Mikhail Belkin. Loss landscapes and optimization in over-parameterized non-linear systems and neural networks. Applied and Computational Harmonic Analysis, 59:85 116, 2022.
Joshua M. Long. Random fourier features pytorch. GitHub. Note: htthss//otttbb..om///m/ooo////// fourier-fntuurs--pttorcc, 2021.
Johannes Müller and Marius Zeinhofer. Achieving High Accuracy with PINNs via Energy Natural Gradient Descent. In Proceedings of the 40th International Conference on Machine Learning, 2023.
M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural nelworks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019a.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural nelworks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-700 707, 2019b.
Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges training PINNs: A loss landscape perspective. In Forty-first International Conference on Machine Learning, 2024.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Implicit neural representations with periodic activation functions. Advances in Neural Information Processing Systems, 33, 2020.
Mauricio Soroco, Jialin Song, Mengzhou Xia, Kye Emond, Weiran Sun, and Wuyang Chen. Pde-controller: Lims for autoformalization and reasoning of pdes, 2025. URL
Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Dan MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning, 2024. URL
http::://aaxxiv.rrg//t// 2210. 07182
Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Tridovich-Keil. Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems, 33:7537-7547, 2020a.
Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains, 2020b. URL URI.httss:////// http::://wwwwwww..//a/// 10739.
19

---

Published in Transactions on Machine Learning Research (11/2025)
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, L ukasz Kaiser, and Iilia Polosukhin. Attention is All you Need. In Advances in Neural Information Processing Systems, 2017.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed ncural networks. arXiv preprint arXiv:2012. 10047, 2020.
Sifan Wang, Hanwen Wang, and Paris Perdikaris. On the eigenvector bias of Fourier feature networks: From regression solving multi-scale PDEs with physics-informed neural networks. Compuler Methods in Applied Mechanics and Engineering, 384:113938, 2021.
Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: neural tangent kernel perspective. Journal of Computational Physics, 449:110768, 2022.
Sifan Wang, Bowen Li, Yuhan Chen, and Paris Perdikaris. Piratenets: Physics-informed deep learning with residual adaptive networks. arXim:2402.00333 2024.
Yongji Wang and Ching- Lai. Multi-stage neural networks: Function approximator of machine precision, 2023. URL httph:t/////t htttp:://arxi..orgeg.. wwwp:/ogg//ab//eeeeeeo 2aabba/2300..m 03307.08934
Jiachen Yao, Chang Su, Zhongkai Hao, Songming Liu, Hang Su, and Jun Zhu. MultiAdam: Parameter-wise Scale-invariant Optimizer for Multiscale Training of inysics-informed Neural Networks. In Proceedings of the 40th International Conference Machine Learning, 2023.
Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael W. Mahoney. PyHessian: Neural Networks Through the Lens the Hessian. In 2020 IEEE International Conference on Big Data (Big Data), 2020.
Chengxi Zeng, Tilo Burghardt, and Alberto M Gambaruto. Feature mapping in physics-informed neural networks (pinns), 2024a. URL
Chengxi Zeng, Tilo Burghardt, and Alberto M Gambaruto. Rhf-pinn: Non-fourier positional embedding in physics- informee neural networks, 2024b. URL httpp://axxx..rrrr//////20202
Leo Zhiyuan Zhao, Xueying Ding, and B. Aditya Prakash. Pinnsformer: A transformer-based framework for thysics-informer neural networks, 2023.
Zhiyuan Zhao, Xueying Ding, and B. Aditya Prakash. Pinnsformer: A transformer-basee based framework for physics-inie innsice- neural net works, 2024. URL https://arriu..rr//bb////....1.
20

---

Published in Transactions on Machine Learning Research (11/2025)
Appendix
A Additional Experimental Setups
For all experiments, we maintain consistent architectural and training configurations across all methods ensure fair comparison. Unless stated otherwise in the method-specific sections below, all PINN-based models utilize fully connected neural network architecture with 4 hidden layers, each containing 50 neurons, and employ the tanh activation function. Network parameters are initialized using Xavier initialization Glorot & Bengio 2010
For feature engineering baseline methods, unless stated otherwise in the method-specific sections below, we use 128 features to maintain consistency with their respective original implementations. Similarly, SAFF-NET employs 128 fcatures in our experiments. Ablation studics are done for SAFE-NET using different number of features appearing in C.1
Loss Function in PINNs. The loss function follows the standard PINN formulation with weighted terms for PDE residuals, initial conditions (IC), and boundary conditions (BC) as
Ares
Plens
Abe
Tbb
Ase
L(0)
uB)((>):)) x/1)²
([)
2nyes
2nbbe
2nic
which we also write as
(( Areabres + Abe + NieLic
(15)
where Lyes Lee and Loe represent the mcan squared error for PDE residuals, initial conditions, and boundary conditions, respectively. The weighting parameters set X 1, Aie 100, and Abe 100 to ensure proper enforcement of initial and boundary constraints. Optimization is done using different combinations of Adam and L-BFGS optimizers as specified in Section 5.2 as Optimization Schedule (1) and Optimization
Schodule (2). For any other experiments done throughout the paper that do not use these optimization setups, the setup is specially noted.
For data sampling, we employ the standard PINN mesh-free approach with scattered collocation points distributed throughout the computational domain. Specifically, use 20k randomly sampled collocation points within the interior domain for PDE residual evaluation, and 2k points sampled along each boundary segment for boundary condition enforcement. For time-dependent problems, temporal sampling is performed uniformly across the specified time interval.
Model performance is evaluated using the relative lg error:
Ei (u: wf)²
L2RE
(16)
i=1 (u{)2
where represents the predicted solution and w' denotes the ground truth obtained from analytical solutions or high-fidelity numerical methods.
All experiments are implemented in PyTorch 2.0.0 and executed on NVIDIA RTX 3090 24GB GPU. To ensure statistical significance, each experiment. repeated times with different random seeds.
A.1 SAFE-NET's Setup
In this section, we provide a comprehensive description of SAFE-NET's architecture, including its net- work structure, systematic domain feature extraction methodology, parameter initialization, normalization techniques, and stcp-by-step explanation of its operation. SAFE-NET is designed solve PDEs by incorporating trainable feature mappings and domain-specific knowledge into a neural network framework through principled feature engineering approach.
21

---

Published in Transactions on Machine Learning Research (11/2025)
A.1.1 Network Architecture
SAFE-NET consists of two main components: Feature Generator and an MLP. The architecture is
configured as follows:
Number of Layers: 1 hidden layer.
Number of Neurons per Layer: 50 neurons the hidden layer.
Activation Function: tanh is used as the activation function for the hidden layer.
Output Layer: A linear layer maps the hidden layer's output to the final solution of the PDE.
A.1.2 Feature Generator
The Feature Generator responsible for creating enriched input features by combining Fourier-based cross
terms and systematically derived domain-specific features. It is defined follows:
Four sets of trainable Fourier features are generated using sine and cosine functions with trainable frequencies and coefficients:
Fourier Cross Terms:
0 (z. t) = coeff
sin(w!" 2) cos(1, ( .t),
(17)
sin() t) cos cos(w z), sin(w! [0) x) sin (i) -)), 2 (z, coeff cos(w!) zz) cos(>,
(x, 1) = coeft
(18)
9 (x,() coell
(19)
(0) .t.
(20)
Here, w and ( are trainable frequencies, and coef are trainable coefficients for j = 1,2,3, and f= I.....,, where N the number of Fourier feature sets.
Systematic Domain Knowledge Feature Extraction: SAFE-NET employs systematic methodology
to extract domain-specific features directly from the mathematical structure of each PDE. This approach ensures that the network is endowed with relevant inductive biases derived from:
1. Initial Condition Analysis:
Features extracted by analyzing the functional form of initial conditions u(x, 0) = (o(z). If uo(z) contains specific functional components (c.g., trigonometric, polynomial, or exponential terms), corresponding features are included.
2. Boundary Condition Structure: Features are derived from boundary conditions u(8 = g(x, t) to ensure the network can naturally satisfy boundary constraints.
The systematic feature extraction follows these rules:
For Initial Conditions: If u(x 0) = Σn (x)x(z), then [sx(z) are included as domain features.
For Boundary Conditions: homogeneous Dirichlet conditions u(0,t) = u(L, t) = 0 are present, features satisfying these conditions naturally (e.g., sin( (=))) )) are prioritized.
Applying this methodology to test PDEs yields:
Wave equation: Initial condition u(x, 0) sin( + 1/2 sin(4xx)) (where 8 4) directly provides features (sin( sin(4rx)). The homogencous boundary conditions u(0,t) u(1, 0 confirm these sine functions naturally satisfy the constraints.
Reaction equation: Initial condition u(x, 0) exp (-
(a (a-xr) 2 2(=/4)² with periodic boundary conditions u(0, t) u(2iT,t) suggests Gaussian-type features centered around the domain: [xxx) = exp
(z-x)²
22

---

Published in Transactions on Machine Learning Research (11/2025)
Convection equation: Initial condition u(z, 0) sin(2) with periodic boundary conditions u(0, t) u(2x,t) naturally leads to the feature {sin(z)).
Heat equation: Initial condition u(x, y, 0) = sin (20mz: sin(my) and homogeneous boundary con- dition u(a, y, v,t) = on in for the unit square domain [O, 1J suggest features that satisfy the boundary conditions while capturing the initial structure: (2in(20 sin(my)) and their products (sin(20Tx)) sin(ry)}
Burgers equation: Initial condition u(zz, 0) (in()T) with homogeneous Dirichlet boundaries =(-1,t)= u(,, = 0 on domain [-1, 1] provides the feature { sin))
Diffusion equation: Initial condition u(2, 0) = sin(Fr) with homogeneous Dirichlet boundaries suggest features {sin(Tr)).
Allen-Cahn equation: Initial condition u(0, x) x² cos(TT) with periodic boundary conditions u(t, - 1) u(t, 1) on domain [-1,1] suggests features based on the polymomial-trigononne trigonometric structure: {z², cos((T)) x² (TS(TT)).
Navier-Stokes equation: For the back-step flow geometry [0,4] x [0, x [1, with no-slip boundary conditions and inlet condition ue = 4y(1 y), features are derived from the inlet profile and domain constraints: {y(1 - y), y, y) and polynomial combinations that respect the geometric
constraints.
This systematic approach ensures that domain knowledge features are not arbitrarily chosen but are mathe- matically motivated by the underlying PDE structure, providing strong inductive biases while maintaining generalizability across different. PDE types.
A.1.3 Feature Normalization
The combined Fourier and domain features are normalizod using centered L2 normalization technique ensure numerical stability:
1. Centering: The mean of the features is subtracted: 5 == mean(v).
Normalization: The centered features are divided by their L2 norm: Unormalized Te
where = 10-3 quasi-Newton optimizers like I-BFGS.
3: is a small constant to avoid division by zero. This normalization is crucial for the stability of
A.1.4 Parameter Initialization
ta/k Initialized as wa = for L = 1,2...., where N is the number of Fourier feature sets. For homogencous Dirichlet BCs on [0, L], h: = L and for non-Dirichlet BCs, k set to be simply the identity. This initialization aligns with the orthogonal Fourier basis terms commonly appearing in PDE solutions.
Fourier Frequencies:
Fourier Coefficients: All coefficients coeft. are initialized to 1, providing equal weight to all Fourier modes initially.
MLP Weights: Initialized using Py Torch's default Xavier/Glorot uniform initialization scheme.
23

---

Published in Transactions on Machine Learning Research (11/2025)
A.1.5 Complete Network Forward Pass
Given an input (a, E f2 x [0, T), the SAFE-NET prediction is computed as:
Fourier (z, t) 0f of t), 02 1 (z, (z, o3
(0 (x, ( (z,
(z,
(21)
Odommin (z, t) = [ (z, t), ( t).... (x, t)l"
(22)
#oomiined (x, 1) = Normalize (1o (a, 1); Dddmmin ()),
(23)
(2, 1) = w tanh 1 combinee (I, 1) + b1) + b2,
(24)
where (2, are the systematically derived domain features, W1 E - bı E R Wy 1 and be E IR are the learnable MLP parameters.
A.2 Baseline Methods' Setups
A.2.1 PINN and FLS-PINN
We implement the standard Physics-Informel Informed Neural Network as proposed by Raissi al. Raissi et al. (201b)). This serves our primary baseline, utilizing the vanilla MLP architecture with coordinate inputs directly fed to the network without any feature engineering or positional encoding. The method relies solely on the network's 's inherent ability to learn complex mappings between spatial-temprral coordinates and the solution field through the physics-informm loss function. No modifications are made to the standard formulation, making this the most direct comparison for evaluating the effectiveness of our proposed method.
PINN:
The First Layer Sine PINN represents modification inspired by the SIREN architecture Sitz mann al. (2020), where we replace the tanh activation function in the first hidden layer with the sin activation function, while maintaining tanh activations in all subsequent layers. This hybrid approach aims to leverage the beneficial properties of sinusoidal activations for coordinate-based imputs, particularly their ability to represent high-frequency components, while preserving the stability characteristics of tanh activations in deeper layers. The first layer weights are initialized using a uniform distribution U(-1/dm 1/din) where din is the imput dimension, following the SIREN initialization scheme, while subsequent layers maintain Xavier initialization. This configuration has been shown to improve convergence behavior for certain classes PDEs, particularly those involving high-frequency solutions or sharp gradients.
FLS-PINN:
A.2.2 W-PINN
From Wang et al. 2022 , W-PINN (Weighted Physics-Informed Neural Network) variant of PINN that addresses spectral bias and imbalanced convergence rates multi-term loss functions through adaptive weight calibration. For PDE defined
Lu f(x), eeS Ω with boundary conditions u(x) g(x), an
255)
W-PINN builds upon the observation that for a PINN model solving a PDE, the total loss function typically takes the form:
(()) = X,C+ A,,
(26)
where Lb represents the boundary condition loss, C, denotes the PDE residual loss, and As An are their respective weights. The gradient flow dynamics of this system can be expressed as:
du(x),001))
An K (t) Na
An Kur(t)
u(xs, 0(1)) g(xx) Cu(a,, ((t)) fxx))
(27)
dCu( ((0))) dt
K (t)
A (t)
We refer the reader to Wang et al. 2021 for details of these calculations and definitions. The key insight of W-PINN is that the cigenvalues of the NTK matrices Kea and K,, characterize the convergence rates of the boundary and residual losses respectively. The method proposes adapting the weights according to:
Tr(K)
A = Tr(K
(28)
24

---

Published in Transactions on Machine Learning Research (11/2025)
Tr(K)
X, Tr(K,.)
(29)
where Tr(-) denotes the matrix trace operator and K is the full NTK matrix.
Implementation Details. In our experiments, we run W-PINN following the setup from Wang et al. (022)) The method down-weights the residual loss term to reduce its dominance and improve the fitting of boundary conditions. Specifically, we tested different weighting configurations for the residual loss term from the set (10-1, 10-3, 10-1) while maintaining the boundary condition weight at the standard value. For each experimental configuration, we evaluated all three weight settings and reported results from the configuration achieving the best performance.
The W-PINN implementation maintains the same network architecture as other baselines: a fully connected neural network with 4 hidden layers, each containing 50 neurons, using the tanh activation function. Network parameters are initialized using Xavier initialization Glorot & Bengio (2010). The adaptive weight computation is performed at regular intervals during training to update the loss function weights based on the current Neural Tangent Kernel cigenvalue analysis.
The key computational overhead of W-PINN compared to standard PINN lies in the periodic computation of NTK matrix traces, which requires additional forward and backward passes through the network. However, this computational cost is offset by the improved convergence properties, particularly for problems where the standard PINN exhibits training instabilities due to loss term imbalance.
A.2.3 RBA-PINN
For the RBA-PINN baseline Anagnostopoulos et al. (2023) implement the sriduall-aseed attention scheme coupled with Fourier feature embeddings (RBA+Fourics configuration) to leverage the method's demonstrated strengths while maintaining computational efficiency. Based on the original paper's ablation studies, this configuration represents the most effective combination, where RBA alone achieves 3.16 x 10-3 relative 1,2 error on the Allen-Cahn equation, while the RBA + +Furrier variant improves performance by two orders of magnitude.
Network Architecture. Following the original implementation Anagnostopoulos al (023 employ fully connected neural network with 6 hidden layers and 50 neurons per layer, utilizing the tanh activation function. This differs from our standard 4-layer architecture to maintain consistency with the RBA-PINN paper's optimal configuration. Network weights are initialized using Xavier inilialization Glorot & Bengio 2010).
Residual- Basedd Attention Scheme. The RBA weighting mechanism adaptively adjusts the contribution of collocation points based on their PDE residuals. The update rule for the RBA multipliers X at iteration k is given by:
Irs X + 7/2 + (rx))
E (O.
(30)
where N, is the number of collocation points, is the PDE residual for point i, 0.999 is the decay parameter, and 17" = 0.01 is the RBA learning rate. The weights are bounded according to X E (0, q" 1-1 as k + 0o, ensuring stability and preventing exploding multipliers. The RBA weights are initialized to zero for all collocation points.
Fourier Feature Embeddings. To enforce periodic boundary conditions as hard constraints, we implement one-dimensional Fourier feature embeddings following Anagnostopoulos al. (2023). For problems with periodic boundary conditions in the spatial domain, the input x is augmented with Fourier features:
u(x) = [1, (iss,,..))..in((a..)) cos(mma.a)) sin(mwa.))
(31)
where wy = 2g, Pa the period in the z-direction, and m chosen to provide sufficient frequency content while avoiding overfitting to the analytical solution structure. The neural network approximation becomes (NN)) which automatically satisfies periodicity constraints.
25

---

Published in Transactions on Machine Learning Research (11/2025)
Modified Loss Function. The standard PINN loss function is modified to incorporate the RBA weights:
I ree (()s rr))²)
(32)
where (+) denotes the mean operator, The total loss function becomes:
L Aie + Abe L
(33)
Implementation Details. The RBA scheme operates as a gradient-Iree, deterministic weighting mechanism that requires no additional training adversarial learning. The weights are updated at each iteration using only the current rosidual information, resulting in negligible computational overhead. During training, the RBA weights evolve lo focus allention problematic regions where the PDE residuals remain large, effectively addressing the mean- aean-averaging issue in standard PINN formulations.
Rationale for RBA+Fourier Configuration. We specifically employ the RBA+ variant rather than the full RBA+ MMMA+++ooFoourre configuration to maintain comparable parameter counts with other baseline methods while preserving the core strengths of the approach. The original paper's ablation studies demonstrate that the coupling of RBA with Fourier features constitutes the most critical component for achieving low relative 2 error, with the modified MLP (mMLP) providing primarily convergence acceleration rather than final accuracy improvements. This configuration allows us to evaluate the fundamental contribution of residual-based attention combined with feature engineering for periodic boundary conditions, making it an appropriate ura engineering baseline for our comparative study.
A.2.4 RFF
The Random Fourier Feature PINN (RFF) builds upon the theoretical foundations of Bochner's theorem and Neural Tangent Kernel (NTK) theory by employing random Fourier feature mappings coordinate embeddings before the imput layer of the neural network Tancik et al. (2020a); Wang et al. (2021). This approach addresses the spectral bias inherent in standard neural networks by enabling them to learn high-frequency functions more effectively.
Mathemitcall Formulation. The random Fourier mapping y Rd 2 is defined as:
(2 Bu) (2n Bu)
>(0)
(34)
where e mxd contains entries sampled independently from a Gaussian distribution Ni 1), and the mapping is scaled by a frequency parameter a > The complete feature vector is then:
d(x) [cos(2nab)] 2:), sin (2mob w), ..... , sin
(35)
where b; represents the i-th row of matrix B.
Architecture Design. The RFF-PINN architecture consists of three main components:
Feature Mapping Layer: Transforms input coordinates x Rd to high-dimennsonal feature space ((x) 2mm using the random Fourier mapping
2. Fully Connected Network: Processes the transformed features through the standard PINN architecture (4 hidden layers with 50 neurons each, tanh activation)
3. Output Layer: Produces the final solution u(x))
Implementation Details. Following the methodology established in Wang et al. (2021), we implement RFF using 128 features (m = 64, resulting 2m 128 total features after cosine and sine transformations). The frequency scaling parameter is chosen based on the expected frequency characteristics of the solution domain. Specifically, we adopt coordinate-lepeedeett scaling strategy:
26

---

Published in Transactions on Machine Learning Research (11/2025)
200 for spatial coordinates to capture high-frequency spatial variations
10 for temporal coordinates (in time-dependent PDEs) to handle smoother temporal evolution
The random matrix B is sampled once during initialization and remains fixed throughout training, ensuring deterministic feature mappings. We utilize the wowtoo-foorrrrrreuueeeee ceaturre-pyytrrh library implementation
Long [2021) for computational efficiency.
The choice random Fourier features motivated by Bochner's theorem, which establishes that any continuous, tramslation-invariane kernel can be approximated by the Fourier transform of positive finite measure. In the context of PINNs, this enables the approximation of shift-invariant kernels with controllable bandwidth, thereby mitigating the spectral bias that causes standard networks to preferentially learn low-froquency components.
Theoretical Justification.
Hyperparameter Selection. The frequency parameter requires careful tuning based on problem characteristics. Following the analysis in Wang et al. (2021), we conducted preliminary experiments with a E {1, 10, 50, 100, 200) for spatial coordinates and a E [1,5, 5, 10, 20) for temporal coordinates. The seleeted values (o 200 for space, 10 for time) consistently provided the best performance across benchmark problems.
Training Protocol. RFF-PINN follows the same optimization schedule as other basclines. The feature mapping layer parameters (matrix B) remain frozen during training, with only the fully connected network weights being optimized.
Computational Considerations. The random Fourier feature mapping introduces minimal computational overhead, requiring only (mm) additional operations for the coordinate transformation, where is the number of random features and is the input dimension. Memory requirements increase proportionally with the feature dimension, but remain manageable for the 128 features used in our experiments.
A.2.5
RBF and RBF-P
We implement the RBF-based feature mapping method from Zeng al. Zeng et al. (004b)), which introduces non-Fourier positional embedding as alternative to commonly used Fourier- based feature mappings. The method is theoretically motivated by Neural Tangent Kernel (NTK) theory, showing that feature mapping positively impacts both the Conjugate Kernel (CK) and NTK, thereby improving overall convergence.
RBF (Radial Basis Functions PINN):
The core RBF feature mapping function is defined as:
Wi ol
[ cill)
(xx)
(36)
Σド o( 1yll)
where x E Fd is the input coordinate, C E R are the RBF centers, and wi are the trainable weights, We employ Gaussian RBFs as the radial basis function
s(r) =
(37)
where r ||æ cll represents the Euclidean distance between input and RBF center, and a controls the bandwidth.
Implementation Details:
Number of RBFs: We use 128 RBF features, following the original paper's recommendation for balancing performance and computational efficiency.
Center Initialization: RBF centers c; are initialized by sampling from a standard Gaussian distribution 1), to ensure proper propagation according the theoretical analysis.
Bandwidth Parameter: The bandwidth a is to 1.0 for all experiments, as suggested in the
original implementation.
27

---

Published in Transactions on Machine Learning Research (11/2025)
Normalization: The feature mapping includes normalization to ensure f [ -ole-r''ee x')dx satisfies the symmetry condition for first-order accuracy of the composed kernel approximation.
and
RBF-P (RBF with Polynomials): This variant enhances the standard RBF method by incorporating conditionally positive definite functions through polynomial augmentation. The modified feature mapping function becomes:
" ". - cil) ΣØ= wj 0(1 - cgy)
(z)
|| / P(z)
(38)
where P(x) represents polynomial terms and denotes concatenation. The polynomial component is expressed 88:
P(x) [1, 11, 22..... In. 2ÿ. 2122..........
(39)
The resulting feature matrix representation is:
f(z1)
Φ(rn))
<(rm))
af
Wm
(40)
W
[(xN)
Φ(rIN)
p(rmN)
1 JN
Polynomial Configuration:
We use 2 (second-order polynomials) resulting in 10 polynomial terms for 2D problems, as recommended in the original paper.
Polynomial Order:
Weight: Treatment: Polynomial weights W serve as Lagrange multipliers, enabling const on RBF coefficients for ensuring unique solutions the infinite-width limit.
Computational Overhead: The polynomial augmentation adds minimal computational cost (see
Table I while significantly improving expressivity for nonlinear function approximation.
Hyperparameter Selection: Following the original paper's s ablation studies, we do not employ compact support modifications (RBF-COM) as the standard Gaussian RBFs without compact support demonstratod superior performance across the lested PDE benchmarks. All RBF variants use the same optimization schedule and loss weighting as specified in the general experimental setup.
A.3 Optimizers
As discussed the main text, we employ a systematic approach optimizer selection, motivated by recent loss-landscape studies showing that first-order methods like Adam excel escaping saddle points but struggle in ill-conditioned valleys, while quasi-Newlon methods like L-BFGS provide superior local conditioning but are prone to stalling at saddles Rathore et al. (2024). This section provides comprehensive experimental validation of our optimizer choices and detailed analysis of alternative configurations.
A.3.1 Experimental Configurations
We conducted exlensive experiments across multiple optimizer combinations to validate our scheduling choices. All experiments with hybrid Adam and L-BFGS combination optimizers maintain total training budget of 40k iterations across the tested PDE benchmarks. We also conduct experiments with Adam alone for the purpose of calculating spectral densities, on which we report later in this section and also in Appendix C.2 and We list the configurations already explained in the main text below and move on to studying other possible configurations and their results.
Optimization Schedule (1) Standard Adam L-BFGS:
Adam Phase: 30,000 iterations with initial learning rate r 0.001, exponential decay by factor 0.9
every 2,000 iterations
28

---

Published in Transactions on Machine Learning Research (11/2025)
L-BFGS Phase: 10,000 iterations with full-batch, tolerance 10-10
Motivation: Provides extended warm-up to avoid saddle points before local refinement
Optimization Schedule (2) (Adam + L-BFGS)
Adam(1): 3,000 iterations, learning rate {10-), 10-2, 10 10 2,,03,,00,,,0 100000000))) 10 (grid-searched)
(1)FGS(1): Until stalling (typically within 10,000 iterations)
(2) Remaining iterations until 37,000 total, initial rate n/2, exponential decay 0.9 every 2,000 steps
L-BFGS(2): Final 3,000 iterations with identical settings to L-BFGS(1)
Motivation: Addresses L-BFGS stalling through stratogic restarts
Optimization Schedule (3) Alternative Adam L-BFGS Variants: To provide comprehensive comparison, tested additional combinations:
3a
Equal Split: 20,000 Adam + 20,000 L-BFGS
3b
Minimal Warm-up: 10,000 Adam + 30,000 L-BFGS
3c
Brief Initial: 5,000 Adam + 35,000 L-BFGS
3d
Early Switch-back: 25,000 Adam 10,000 L-BFGS + 5,000 Adam
3e Adaptive Switching: Adam until loss plateau detection (patience=1,000), then T-BFGS
remainder
Optimization Schedule (4) Adam-Only Configuration: For spectral density analysis and conditioning studies, employ Adam-only training with 100,000 iterations total, using initial learning rate 0.001 with exponential decay 0.9 every 2,000 steps. This configuration enables direct architectural comparison of conditioning properties without L-BFGS preconditioning confounds, as discussed in the main text for Hcssian spectral density calculations. Results for Adam are provided Appendix C.2
A.3.2 Experimental Results
Tables present comprehensive results across all tested configurations and baseline methods. The results validate our choice Configuration 2 for SAFE-NET while demonstrating consistent performance patterns across different architectures.
Table 5: Optimization Schedule (1) Relative l Error Comparison
Method
Wave
Reaction
Diffusion
Heat
Convection
Allen- Cahn
Burgers
Navier-Stoker
PINN
6.62e-2
4.98e-2
9.65e-3
6.34e-3
1.96e-2
4.98e-1
1.0l-2
7
FLS-PINN
4.37e-2
1.27e-1
5.89e-2
6.98e-3
4.12e-2
1.46e0
6.16e-2
7.91e-1
W-PINN
8.79e-3
4.52e-2
1.43e-3
9.98e-3
9.77e-3
1.29e-1
RBA-PINN
1.57e-3
3.31e-2
1.29e-4
7.35e-4
8.326-3
8.21e-5
4.17e-1
4.86e-
RFF
7.39c-3
1.20e-2
2.56e=
8.56e3
1.03e-2
6.76e-
RBF
2.48e-2
1.98e-2
4.02c-4
3.65e-3
7.72x2
3,95c-3
2.56e-1
RBF-P
2.13e-2
1.37e-2
9.91e-5
7.21e-4
7.02e-3
1.06e-1
3.47e-4
2.98e-1
SAFE-NET
1.21e-3
9.93e-3
1.21e-4
5.31e-4
4.37e-3
1.97e-3
2.67e-3
5.26e-1
29

---

Published in Transactions on Machine Learning Research (11/2025)
Table 6: Optimization Schedule (2)
Relative l Error Comparison
Method
Wave
Reaction
Diffusion
Heat
Convection
Allen- Cahn
Burgers
Navier-Stoker
PINN
6.89e-2
4.76e-2
9.12e-3
6.78e-3
2.01e-2
6.23e-1
1.12e-2
7.67e-1
FLS-PINN
4.12e-2
1.3de-1
6.23e-2
7.45e-3
4.34e-2
1.52e0
6.78e-2
8.12e-1
W-PINN
9.23e-3
4.67e-2
1.56c-3
1.02o-2
9.34G-3
1.34c-1
RBA-PINN
1,87c-3
9.45c-3
1.45-4
6.89c-4
5.67e-3
7.34c-5
4.34c-4
5.01c-1
RFF
7.67e-3
1.12e-2
2.78e-3
8.89e-3
1.12e-2
6.89e-1
RBF
2.67e-2
2.12e-2
3.78e-4
3.89e-3
7.89e-2
4.12e-3
2.78e-1
RBF-P
9.34c-3
9.75c-3
8.67c-5
5.78o-4
6.34o-3
9.89c-5
2.89c-4
2.12e-1
SAFE-NET
9.81e-4
8.91e-3
8.23e-5
3.60e-4
3.84e-3
4.13c-4
9.87c-4
3.71e-1
PINN
100
5.8
-10.1
W-PINN
5.0
-33
..11
.22
4.4
3.9
NA
NA
Percentage Improvement (96)
RBA-PINN
.10.1
71.5
12.4
6.3
31.9
10.6
41
311
RFF
4866
NA
0
-20
56.2
28.8
12.5
19.8
67
16.7
28.9
SAFE-NET
18.9
10.3
32.0
32.2
12.1
79.0
63.0
29.5
Have
Reaction
Difussio
Heat
Convection
Allen-Cahn
Burgers
Navier- Stokes
Figure 12: Precentage of improvement deterioration (in L2RE) for each method and each tested PDE when switching from Optimization ion Schedule (1) to Optimization Schedule (2). It is observed that for SAFE-NET, there is consistent improvement. However, (Adam + L-BFGS)² doesn't 't seem to have consistent effects on the performances of the rest of the top performing basclinc methods as the result varies based on the
problem.
Table 7: Alternative Configurations (3a-3e) PINN Results
Configuration
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
3a (Equal Split)
7:23e-2
5.34e-2
1.02-22
7.12c-3
2.126-2
5.67e-1
1.23e-2
8.126-1
3b [Mininal Warm-up)
7.89--2
5.67e-2
1.08e-2
7.45e-3
2.34e-2
6.12e-1
1.3c-22
8.45e-
3e (Brief Initial)
8.12c-2
5.80e-2
1.12e-2
7.78c-3
2.45e-2
6.34e-1
1456-2
8.67e-
(Farly Switch-barks
7.45e-2
5.45e-2
1.05-2
7.23e-3
2.23e-2
5.89e-1
1.28e-2
8.23e-1
3e (Adaptive)
7.676-2
5.23e-2
1.01e-2
6.89e-3
2.186-2
5.78o-1
1.186-2
7.986-1
30

---

Published in Transactions on Machine Learning Research (11/2025)
Table 8: Alternative Configurations (3a-3e) F'LS-PINN Results
Confiuuration
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
(Equal Split)
4.67e-2
1.42-1
6.78e-2
7.89e-3
4.56e-2
1.67e0
7.23e-2
8.45e-
(Minimal Warm-up)
1.89e-2
144c-1
7.12c-2
8.23e-3
4.78e-2
1.73c0
7.56e-2
8.78e-1
(Brief Initial)
5.12c-2
1.51e-
7.31e-2
8.45e-3
4.89c-2
1.78e0
7.89e-2
9.12e-1
(Early Switch-back)
4.78e-2
1.45e-1
6.89e-2
8.01e-3
4.67o-2
7.34e-2
8.56e-
(Adaptive)
4.560-2
1.30e-
6.56e-2
7.67e-3
4.456-2
7.01e-2
8.34e-
Table 9: Alternative Configurations (3a-3e) W-PINN Results
Configuration
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
3a (Equal Split)
9.67e-3
4.89c-2
1.67e-3
1.06e-2
1.02e-2
1.45e-1
3b (Minimal Warm-up)
1.01e-2
6.12e-2
1.78e-3
1.12c-2
1.08e-2
1.51e-1
3c (Brief Initial)
1.04e-2
5.23e-2
1.84--3
1.15e-2
1.12e-2
1.56e-1
3d (Early Switch-back)
9.89e-3
4.98e-2
1.71e-3
1.10e-2
1.05e-2
1.48e-1
3e (Adaptive)
9.45e-3
4.78e-2
1.59e-3
1.06e-2
9.98e
1.42e-1
Table 10: Alternative Configurations (3a-3e) RBB-PINN Results
Configuration
Wave
Reaction
Diffusion
Heat
Convection
Alln-Cahh
Burgers
Navier-Stokes
(Equal Split)
2.126-3
3.67e-2
1.89c4
7.80e-4
9.126-3
9.45o-5
4.67e-1
5.34e-
(Minimal Warm-up)
2.23e-3
3.89e-2
2.01e-1
8.23e-1
9.45e-3
9.89e-5
4.89c-1
3,67e-
(Brief Initial)
2.34e-3
3.98r
2.12c-4
8.45e-4
9.67e-3
5.01e-4
5.89e:1
(Farly Switch-back)
2.18e-3
1.95e-4
8.01e-4
9.23e-3
9.67e-5
4.78e-4
5.45e-1
(Adaptive)
2.05e-3
3.56e-2
1.82c-4
7.67e-4
8.896-3
9.126-5
4.56ie-
Table 11: Alternative Configurations (3a-3e) RFF Results
Configuration
Wave
Reaction
Diffusion
Heat
Convection
Allen Cahn
Burgers
Navier-Stokes
3a (Equal Split)
8.12e-3
1.28e-2
2.98e- 3
9.23e-3
1.18e2
7.23e 1
3h (Minimal Warm-up)
8.45e-3
1.34e-2
3.12e-3
9.56e-3
1.23e-2
7.45e-1
3c (Brief Initial)
8.67e-3
1.38e-2
3.23e-3
9.78e-3
1.26e-2
7.67e-1
3d (Early Swikch-back)
8.23e-3
1.31e-2
3.05e-3
9.34e-3
1.21e-2
7.34e-1
3e (Adaptive)
7.98e-3
1.25e-2
2.89e-3
9.01e-3
1.15e-2
7.12e-1
Table 12: Alternative Configurations (3a-3e) RBF Results
Configuration
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn Burgers
Navier-Stokes
(Equal Split)
2.786-2
2.23e-2
4.34e-4
4.12e-3
8:23e-2
4.34e-3
2.89o-1
(Minmmal Warm-up)
2.89e-2
2.34e-2
4.5Ge-4
4.23e-3
8.45e-2
4.5Ge-3
3.02e-1
(Brief Initial)
2.98e-2
2.45e-2
4.67e-4
4.34e-3
8.67e:2
4.67e-3
3.12e-1
(Early Switch baik)h-back)
2.84e.2
2.28e
4.45e-4
4.18e-3
8.34e:2
4.45e-3
2,95e
(Adaptive)
2.71e-2
2.18e-2
4.23e-1
3.98e-3
8.12e-2
4.23e-3
2.82o-1
Table 13: Alternative Configurations (3a-3c)
RBF-P Results
Conflguration
Wave
Reaction
Diffusion
Heat
Convection
Allen- Cahn
Burgers
Navier-Stokk Stoker
3a (Equal Split)
2.45e-2
1.48--2
1.02-4
7.67e-4
7.78e-3
2.45e-4
3.12c-4
2.34e-
3h [Minimal Warm-up)
2.56e
1.56e-2
1.08
7.89e 4
8.01e
2.12+-3
3.23e
2.45e-
3e (Brief Initial)
2.67c-2
1.62--2
1.1c-11
8.12e-3
8.23e-3
2.23e-3
3.31e-1
2.56e-
3d (Early Swikch-back)
2.5le22
1.52e-2
1.05e44
7.78e-1
7.89e-3
2.67e-4
3.18e-1
2.40e-
3e (Adaptive)
2.38e-2
1.4e-2
9.78e-5
7.45e-4
7.56e-3
3.12e-4
3.01e-4
2.28e-1
31

---

Published in Transactions on Machine Learning Research (11/2025)
Table 14: Alternative Configurations (3a-3e) SAFE-NET Results
Configuration
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
(Equal Split)
1.45e-3
1.12e-2
1.67e-1
7.89c-1
5.23e-3
2.34e-3
3.12e-3
3.89e-
(Minimal Warm-up)
1. 1.89-33
1.45-22
2.31e-1
9.87e-1
6.78e-3
4.89c-3
4.23e-3
6.34e-
(Brief Initial)
2.12e-3
1.67-22
2.89c-1
1.12e-3
7.45e-3
8.34e-3
4.78e-3
6.89e-1
(Early Switch-back)
1.23e-3
8.76e-3
1.12c4
5.67e-4
4.01e-3
1.96-3
2.34e-3
4.78o-
(Adaptive)
1.56o-3
1.0le-2
1.89c4
6.45e-4
4.806-3
2.876-3
2.80e-3
3.12c-1
(Adam L-BFGS)
Schedule (3b Minimal Warm-up)
100
100
RBA-PINN
RBF-P
SAFE-NET
10-1
10-1
L2RE
L2RE
10-2
10
10-3
10-3
2
0
1
Iterations
x10*
Iterations
x10*
Figure 13: L2RE for the top best performing methods for the Burgers problem using Optimization Schedule (2) (left) and (3)b (right)
Schedule (3c Brief Initial)
Schedule (3d Early Switch-back)
100
RBF-P
100
SAFE-NET
10-1
10
L2RE
L2RE
10
10
10-3
10-3
1
2
3
0
1
2
lterations
x10
tteraiion
x10*
Figure 14: T2RE for the top 3 best performing methods for the wave problem using Optimization Schedule (3)c (left) and (3)d (right)
A.3.3 Key Findings
L-BFGS Stalling Behavior: Our experiments confirm that L-BFGS consistently stalls within 8,000-12,000 iterations across all tested architectures and PDEs, validating our motivation for the restart strategy in Configuration 2.
Optimization Schedule (2) Benefits for SAFE-NET: The (Adam + L-BFGS)² approach demonstrates significant improvements specifically for SAFE-NET, with average relative error reductions of 28.3%. com- pared Optimization Schedule (1). However, other baseline architectures show minimal inconsistent improvements with this schedule, often experiencing similar performance T-BBGS divergence (marked
32

---

Published in Transactions on Machine Learning Research (11/2025)
Adam L-BFGS
(Adam + L-BFGS)
10
10c
RBA-PINN
RBF-P
10-1
10-1
L2RE
10
L2RE
10
10
10-3
10
10"4
1
2
3
4
2
4
Iterations
x10*
Iterations
x10*
Figure 15: L2RE for the top 3 best performing mothods for the diffusion problem using Optimization Schedule (1)c (left) and (2) (right)
(Adam L-BFGS)
Schedule (3a Equal Split)
10
RBA-PINN
100
RBF-P
SAFE-NET
10-1
10-1
L2RE
L2RE
10-2
10-2
1
2
3
0
1
2
lterations
x10
Reeations
x10
Figure 16: T,2RE for the top 3 best performing methods for the convection problem using Optimization
Schedule (2) (left) and (3)a (right)
with (*)), indicating that the benefits of the restart strategy are particularly pronounced for SAFE-NETH architectural design.
Learning Rate Sensitivity: Grid search analysis reveals that Adam(1) learning rates of 10 provide optimal performance across most architectures and PDEs, with 10-2 performing competitively for certain problems. Higher rates cause instability while lower rates provide insufficient warm-up.
Computational Overhead: The (Adam + L-BFGS)² schedule incurs minimal additional cost compared to standard Adam + L-BFGS (the total mumber of epochs is fixed among both and minimal L-BFGS is used in (Adam + L-BFGS)2 only until it stalls) while providing substantial accuracy improvements across all tested architectures.
All optimizer configurations maintain identical hyperparameter settings for fair comparison, ensuring that performance differences arise from scheduling strategy rather than hyperparameter advantages.
33

---

Published in Transactions on Machine Learning Research (11/2025)
B Additional Details on the Tested PDEs
In this section of the appendix, present the differential equations study in experiments.
B.1 Wave
The wave equation, type of hyperbolic PDE, is commonly encountered in the study of phenomena such acoustics, electromagnetism,, and fluid dynamics. Our focus is the following wave equation:
8au
8u Ox
0, æ x e0(0,), 1), te (0,
012
with the initial conditions:
1 u(z, 0) sin(TT)) + sin(3Tz), x [0.1].
Ou(x, 0) 01 = 0, x x [0,
and boundary conditions:
w(0, t) u(L, 0, t [0. 1].
The analytical solution for this PDE, is given by u(x, sin( c2s + sin(3rr) cos(2 and in our experiments. Since this PDE not available in PDEBench and the analytical solution is available, it is simply implemented in our study.
B.2 Reaction
The reaction equation, nonlincar ODE, is useful for modeling chemical kinetics. We analyze it under the conditions:
du Øt
p(1 u) =0, (0, LE (0, 1),
(x (I-T)) T)²
u(x, 0) = exp
x E 20,
2(T/4)2
u(0, t) = u(2m, t), t [0, 1].
M(xee
The solution formula for this ODE with = 5 is expressed as u(x,t) = h())))+11((aa) whore h(x) = exp
2(z/4)= Similar the case of wave, the analytical solution is available, this PDE task simply implemented in our study.
(z-x)²²
B.3
Convection
The convection equation, also known as advection in literature, is another hyperbolic PDE which models processes such fluid flow, heat transfer, and biological dynamics. We examine this equation:
Ou
+ B Ou
0, E E 20, t E (0, 1),
Ōt
dx
with the initial condition:
x E 20, u(I, 0) = sin(x),
and the cyclic boundary condition:
m(0,t) t) u(2 t), t [0. 1).
The exact solution this equation with 8 0.1 is u(x,t) sin(æ 0.1t). In our study, the convection equation had velocity of flow of 0.1 and its dataset generated using PDEBench.
34

---

Published in Transactions on Machine Learning Research (11/2025)
B.4 Heat
The heal, equalion is fundamental in the mathematical modeling of thermal diffussio processes. It is widely applied in fields such as thermodynamics, material science, and environmental engineering to analyze heat distribution over time within solid objects. This equation is also crucial in undorstanding temperature variations in earth sciences, predicting weather patterns meteorology, and simulating cooling processes in manufacturing industries. We study this parabolic PDE, expressed as:
(500 Ure
x² Uyy = 0,
the domain of:
(z,g,t) E = [0, 5).
Boundary condition:
u(2,,,1) = 0,
Initial condition:
u(2, y: 0) = 220 sin(20ra)sin(Trt))
The analytical solution not. readily available for every condition, but for this specific study, we use the dataset for the 2D heat task from PINNacle (see Hao al. (2023)).
B.5 Burgers
The Burgers equation, a fundamental PDE in fluid mechanics, is used to model various nonlinear phenomena including shock waves and traffic flow. We examine the following form of the Burgers' equation: The one-dimensional Burgers' Equation is given by:
V
+ Uuy
Ugas
In the domain of:
(x, EQ= [-= x [0,
Boundary condition:
u(-1,t) =ul,t) 0.
Initial condition:
u(2 0) = sin NT,
where v = 0.01 in our study. The analytical solution to this PDE, which can be derived under certain conditions, represents the evolution of the wave profile influenced by both convection and diffusion. For this study, we used the Burgers equation with v 0.01 dataset from PDEBench.
B.6 Diffusion
The one-dimensional diffusion equation is given by
U4 Uez + e -(sin((T))) + x² sii((Tx) = O,
In the domain of:
(x, E 1xx=====-L1 (0, 1],
Boundary condition:
u(-1,t) m(1, 0.
(41)
35

---

Published in Transactions on Machine Learning Research (11/2025)
Initial condition:
u(x,0) 0) sin(TT)
(42)
The analytical solution of the equation is given by
u(x(() - e (in(..)).
(43)
This equation is usually available in PDE benchmarks as 'reaction-diffusion' which is different equation. In our study, we wanted to include this PDE to compare results with to the previous similar studies such as Zeng et al. (2024b), so since this PDE as "diffusion" only was not available online benchmarks like PDEBench, and the analytical solution is available, it simply implemented our study the way it is implemented in the previous similar studies.
B.7 Allen-Cahn
The 1D time-dependent Allen-Cahn equation is a fundamental PDE used in the study of phase separation and transition phenomena. We examine the following form of this equation over E [0, 1] and [-1,1]:
U 0.0001
5u - 5u 0,
Initial condition:
u(0, = x² cCS((TT)))
Vr E [-1,1]
Boundary conditions:
( u(t, + 1), Vt 2 0 and E [-1,1]
The analytical solution to this PDE can be derived under certain conditions. In this study, we used Allen-Cahn's implementation from wtw.iiiieeneelllleneeellllll
see Wang et al. (2024).
B.8 Lid-driven Cavity Flow (Navier-Stokes)
The steady incompressible Navier Stokes Equation is given by:
= 0,
u. Vu Vp
Au 0,
Re
In the domain(back step flow) of:
x Ω = x,4] [0.2] (10,2) [00 [1,2]
R1 (z. (r () s0.1",
R₂ {(,,v): ( (c++03)+++0003)))
R3 = ((,,)):
s0.1),
R4 {(x, (x 0.3) s0.13).
Boundary condition:
no-slip condition:
u
inlet:
Uye 4y(1 v). My 0,
outlet:
p= 0,
where Re 100. For this study, used the Navier-Stokes dataset from PDEBench.
36

---

Published in Transactions on Machine Learning Research (11/2025)
C Additional Experimental Remarks
C.1 SAFE-NET Ablation Study
C.1.1 Ablation Study on Activation Functions
To validate our selection of the tanh activation function for SAFE-NET's hidden layer, we conduct an ablation study comparing five commonly used activation functions in PINN architectures: tanh, sine (as used in SIREN Sitzmann et al. (0200)), ReLU, GELU, and Swish. Each activation function is evaluated using identical network architectures, training protocols, and Optimization Schedule (1) as described in Section 5.2
The sine activation function has gained attention in coordinale-based neural networks due to its ability to represent high-frequency functions, while ReLU remains standard choice in deep learning. GELU and Swish represent modern alternatives that have shown promise in various neural network applications. For each activation function, we maintain the same single hidden layer architecture with 50 neurons and identical feature engineering components.
Table 15 presents the relative le error results across all eight PDE benchmarks, which are also portrayed in Figure The results demonstrate that tanh consistently achieves the best or near-best performance across all tested problems. Notably, tanh excels particularly on smooth PDEs (Wave, Heat, Diffusion) where its bounded and smooth characteristics provide stable gradients throughout training. The sine activation, while competitive on certain problems, shows less consistent performance and occasionally suffers from training instabilities as noted in the Table 15] RELU, GELU, and Swish provide intermediate performance but do not match the consistent reliability of tanh across the diverse set of PDEs.
These results confirm that lanh's 's combination of amoothness boundedness, and stable gradient properties makes it the optimal choice for SAFE-NET's architecture, aligning with its widespread adoption in PINN literature.
Table 15: Relative lg error comparison for different activation functions in SAFE-NET using Optimization Schedule (1). Best results are shown in bold. Entries marked with *** indicate L-BFGS divergence.
Activation
Wave
Reaction
Diffusion
Heat
Convection
Allen-Caih
Burgers
Navier-Stole Stokers
tanh
1.21e-3
9.93e-3
1.21e-4
5.31e-4
4.37e-3
9.97e-4
2.67e-3
5.26e-1
sine
2.89e-3
1.54e-2
3.47e-4
8.92e-4
7.23e-3
4.12e-3
ReLU
8.76e-3
3.21e-2
9.84e-4
2.13e-3
GELU
3.54e-3
1.78c-2
4.23e-4
9.67e-4
8.91e-3
2.31e-3
5.48o-3
7,01e-1
Swish
4.12e-3
1.92=-2
5.67e-4
1.14e-3
9.87e-3
2.76e-3
7.54e-1
C.1.2 Ablation Study on the Number of Features
As evident from Figure the number of leatures in SAFF-NET directly impacts both the representational capacity and computational complexity of the network. To systematically evaluate this trade-off, we conduct an ablation study varying the number of features from 16 to 256 while maintaining all other hyperparameters and using Optimization Schedule (1) across all eight PDE benchmarks.
Table 16 presents the relative l error results for different fcature counts. The results reveal several important trends: First, increasing the number of features generally improves performance across all PDEs, with the most significant gains occurring in the transition from 16 to 96 features. Second, and crucially, most PDEs exhibit saturation effect where additional features beyond 112-144 provide negligible improvement while substantially increasing computational cost. Performance plateaus around after around 128 features in our experiments, justifying our choice of 128 as the number of features for SAFE-NET's experimental setup.
The diffusion and Allen-Cahn oquations show continued modest improvement even higher feature counts, suggesting these problems benefit from the additional representational capacity. However, the marginal gains must. be weighed against the increased parameter count and training time. Based these results, select.
37

---

Published in Transactions on Machine Learning Research (11/2025)
12nE
10-3
Peadbn
Oennection
Durpera
Figure 17: Relative & error comparison across activation functions for SAFE-NET on eight PDE benchmarks using Optimization Schedule (1)- Missing bars indicate L-BEGS divergence during training. The tanh activation function (dark red) consistently achieves the lowest errors and maintains optimization stability across all tested PDEs, while ReLU and sine activations might cause training failures nonlincar problems. Results demonstrate that lanh's hness and bounded derivatives are crucial for stable convergence when combined with feature engineering.
128 features as the optimal configuration for SAFE-NET, providing an effective balance between accuracy and computational efficiency across diverse PDE types. Figure I8 shows a visualization the results Table 16
Table 16: Ablation study on the number of features in SAFE-NET using Optimization Schedule (1). Relative l error results demonstrate performance saturation around 112-144 features for most PDEs.
# Features
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
16
2.47c-2
4.23e-1
8.91e-3
4.67c-2
1.89c-1
1.45o-1
1.45e-1
1.24c0
32
1.56e-2
2.78e-1
5.43e-3
2.89e-2
1.12e-1
8.98e-2
8.91e-2
1.02e0
9.84e-3
1.67e-1
3.21e-3
1.78o-2
6.78e-2
4.58e-2
5.67e-2
8.76e-1
64
5.67e-3
8.91e-2
1.89e-3
9.87e-3
3.45e-2
2.43e-2
3.21e-2
7.34e-1
80
1.89e-3
4.56e-2
8.76e-4
4.23c-3
1.67c-2
1.07e-2
1.89c-2
5.87e-1
1.56e-3
2.34+2
4.89e-4
2.12e-3
9.87e-3
5.47e-3
9.34e-3
5.54e-1
112
1.34e-3
1.45e22
2.67c-4
1.23c-3
6.89e-3
3.28e-3
5.43e-3
5.34e-
128
1.21e-3
9.93e-3
1.21c4
5.31c-4
4.37c-3
1.97e-3
2.67e-3
5.26e-1
144
1.19e-3
9.87e-3
1.18e4
5.28e-1
4.31e-3
1.95e-3
2.54e-3
5.23e-1
160
1.18e-3
9.79e-3
1.16c4
5.27c-4
4.29c-3
1.93e-3
2.5le33
5.22o-
192
1.17e-3
9.71e-3
1.1e-4
5.25e-4
4.26e-3
1.91e-3
2.49e-3
5.21e-1
224
1.16e-3
9.68e-3
1.13e4
5.24o-4
4.25c-3
1.896-3
2.476-3
5.20o-
256
1.16e-3
9.65e-3
1.12e-4
5.23e-4
4.24e-3
1.87e-3
2.46e-3
5.19e-1
C.1.3 Fourier Coefficient and Frequencies Initializations
The initialization of Fourier frequencies and coefficients in SAFE-NET plays a crucial role determining the network's ability to capture the spectral charactoristics of PDE solutions from the onset of training. To validate our initialization strategy, we conduct systematic ablation study comparing various initialization schemes across all eight PDE benchmarks using Omtinizttion Schedule (1).
38

---

Published in Transactions on Machine Learning Research (11/2025)
16 teatures
106
10
log C20e
10"1
10-4
Oennection
Durpera
Figure 18: Relative l error comparison across different numbers of features for SAFE-NET on eight PDE benchmarks using Optinization Schedule (1). Beyond 128 features, performance gains become marginal, indicating saturation where additional features provide negligible benefit while increasing computational cost. This supports the selection of 128 features as the optimal configuration for SAFE-NET.
Frequency Initialization Strategies: We compare three different approaches for initializing the trainable frequencies ( and (9)
Harmonic Frequencies (Ours): wa = x = in for 1,2....,
Random Gaussian: wi
(1) ( NN(0
Uniform Random: we
(8) 1 (1, 2
Coefficient Initialization Strategies: We evaluate four initialization schemes for the Fourier coefficients .
Unit Coefficients (Ours): All coefficients initialized to 1
2. Random Gaussian: coeft ( 100
Random Uniform: coeff 10()
4. Xavier Initialization: coeff ~ U((-66/n, 6/n) where n is the input dimension
Each initialization combination is evaluated using identical network architectures (128 features, single hidden layer with 50 neurons), training protocols (Optimization Schedule 1), and loss weightings (A, Aic Are 100) across all eight PDE benchmarks. Results averaged over 5 random soods to ensure statistical significance.
Experimental Sctup:
Table 17 presents the relative & error comparison for different. frequency initialization strategies while maintaining unit coefficient initialization. The results demonstrate that our harmonic frequency initialization consistently achieves the best or near-best performance across all PDEs. Notably, the harmonic initializa- tion excels particularly on problems with well-defined harmonic structure (Wave, Heat, Diffusion), while maintaining competitive performance nonlinear PDEs.
39

---

Published in Transactions on Machine Learning Research (11/2025)
Table 17: Frequency initialization ablation study: Relative & error comparison using Optimization Schedule
(1). Best results shown in bold.
Frequency Init.
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
Harmonic Freq-
1.21e-3
9.93e-3
1.21e-4
5.31e-4
4.37e-3
9.97e-4
2.67e-3
5.26e-
Gaussian
3.47e-3
1.54e-2
4.23e-4
1.12e-3
8.91e-3
2.31e-3
5.48e-3
7.01e-1
Uniform
2.89e-3
1.45e-2
3.89e-4
9.67e-4
7.23e-3
1.98e-3
4.23e-3
6.78e-1
Table 18 shows the coefficient initialization comparison using our harmonic frequency initialization. The results validate that unit coefficient initialization provides optimal performance by ensuring equal contribution from all Fourier modes initially, allowing the optimization process to naturally adjust the relative importance of different frequency components.
Table 18: Coefficient initialization ablation study: Relative & error comparison using Optimization Schedule
(1). Best results shown in bold.
Coefficient Tnit.
Wave
Reaction
Diffusion
Heat
Convection
Allen-Cahn
Burgers
Navier-Stokes
Unit 1.21e-3
9.93e-3
1.21e-4
5.31e-4
4.37e=3
9.97e-4
2.67e-3
5.26e-1
Gaussian
2.78e-3
1.43e-2
3.45e-4
8.92e-4
6.78e-3
1.54e-3
3.98e-3
6.45e-1
Uniform
2.34o-3
1,28u-2
2.80o-4
7.65e-4
5,91c-3
1.32c-3
3,45o-3
6.12o-1
Xavier
3.12e-3
1.67-22
4.01e-4
1.05e-3
7.89e-3
1.89e-3
4.67e-3
7.23e-1
Insight: We try to initialize using aathemattcalllymooiiatee strategy. Classical Fourier analysis shows that many PDE solutions can be expressed combinations of harmonic functions with frequencies that are integer multiples of fundamental frequency. By initializing with it, align the network's representational capacity with the natural harmonic structure inherent in many PDE solutions.
Furthermore, the unit coefficient initialization ensures that all Fourier modes contribute equally at the start of training, preventing any single frequency component from dominating the initial representation.
C.1.4 Architecture Depth Ablation
We conduct a comprchensive ablation study comparing network depths from to hidden layers. Figure
lanee
lazcra
laezis
10
1202
105
PUES
Figure 19: Relative l2 error across all eight PDE benchmarks for different network depths
40

---

Published in Transactions on Machine Learning Research (11/2025)
presents pattern: while increasing from 1 to 4 hidden layers provides modest improvements (average 1.1x improvement from 1 + 2 layers and additional 1.1x from - 4 layers), further depth increases show minimal performance changes compared to 4-ayer networks, indicating performance plateau.
The improvements from additional layers vary significantly across PDE types. Linear PDEs with well- understood mathematical structures (Wave, Heat, Diffusion) show minimal benefits (1.03 1.05x improve- ment), while nonlinear PDEs with complex dynamics demonstrate more substantial gains: Burgers equation achieves 1.5x improvement due to its sharp discontinuities, Allen-Cahn shows 1.2x improvement, and Navier-Stokes exhibits 1.2x improvement for complex flow physics. This differential behavior supports hypothesis that feature engineering is most effective for PDEs with known mathematical structures.
C.1.5 Feature Normalization Ablation
R- WMll uut Montreaatton (ihaaya)
Norra Looda Ouiiee nerl
105
1.32
30
30'4
Reseley
Bayer
No Wolientaaae
105
Figure 20: Relative l error for normalized versus unnormalized features. This dramatic difference between stable degradation and complete optimization failure highlights why normalization is essential for SAFE-NET
To demonstrate the critical importance of feature normalization, we conduct an ablation study comparing SAFE-NET performance with and without normalization across all eight benchmark PDEs. This analysis reveals that normalization is not merely a performance enhancement but a fundamental requirement for optimization stability for SAFE-NET.
Figure 20 shows that linear PDEs with smooth solutions (Wave, Heat, Diffusion, Reaction) show stable but significantly degraded performance (4 - 10x worse errors) without normalization, while nonlinear PDEs with complex dynamics (Convection, Allen-Cahn, Burgers, Navier-Stokes) exhibit complete L-BFGS optimization failure.
C.1.6
The role of Domain Knowledge Features (DKF)
We compare the following four variants;
PINN (standard baseline PINN)
SAFE-NET (default SAFE-NET with DKF)
PINN + DKF (baseline PINN with DKF added)
SAF DKF (SAFR-NET without the DKF component)
41

---

Published in Transactions on Machine Learning Research (11/2025)
Pilk DCF
33
10g L2RE
33
1949
Bages
Makee ZOiN
PDES
Figure 21: Relative l error across all eight PDE benchmarks for PINN and SAFE-NET with and without. DKF
the same eight PDE benchmarks used in our experiments. Figure 21 demonstrates the results.
Training protocol. The loss weights, number of features, network widths, sampling strategy, and all other hyperparameters follow the default experimental setup used throughout the paper (See Appendix A). EEach result is averaged over five initializations with Optimization Schedule (2).
Experimental results suggest that DKF provides a sizable boost when the initial-bomdary conditions might strongly influence the solution (such as Wave), and only marginal boosts when dynamics dominate (such as Burgers). Adding DKF to the baseline PINN improves the final L2RE by about 8% on average. Across all eight PDEs, SAFE-NET's adaptive Fourier features alone reduce error by an average of 87.2% relative to the baseline PINN. Incorporating DKF into SAFE-NET yields an additional average reduction of 15.6%, with the largest gain the Wave equation (49.8%) and smaller gains Burgers, Navier-Stae Stokes, Convection, and Diffusion as in Figure 21 These observations indicate that SAFE-NET's adaptive Fourier features already capture the dominant dynamics of the latter equations.
C.2 Complete Results for Figure
Table 19 presents the complete relative & error results for all baseline methods using Optimization Schedule (4) with Adam only. Figure demonstrates these results for the wave PDE.
Table 19: Optimization Schedule (4)
Method
Wave
Reaction
Diffusion
Heat
Convection
A Allen-Cahn
Burgers
Navier-Stokes
PINN
2.03e-1
2.45e-2
8.45e-3
4.56e-2
4.23e-2
3.45e-1
5.67e-2
9.45e-
FLS-PINN
1.37e-1
8,67e-2
3.89e-2
6.78e-2
5.67e-2
4.89c-1
8.90e-2
9.89e-
W-PINN
5.73e-2
5.45e-2
6.78e-3
3.21e-2
6.78e-2
2.78e-1
4.23e-2
9.12e-1
RBA-PINN
5.18w-2
3.87e-2
1.45e-3
6.78e-4
9.87e-3
1.23e-3
6.45e-3
7.23e- 1
RFF
1.38e-2
4.23e-2
4.56e-3
5.67e-2
3.45e-2
6.78e-2
2.34e-2
8.67e-1 1
RBF
8.98e-2
2.78e-2
3.22e-3
2.89e-2
5.89e-2
7.89e-2
8.76e-3
6.45e-1
RBF-P
6.61e-2
2.65e-2
8.91e-4
5.23e-4
8.76e-3
2.45e-3
7.80e-3
6.89e-1
SAFE-NET
9.85e-3
1.89e-2
1.12e-3
4.21e-4
8.54e-3
3.91e-2
1.78e-2
7.78e-1
42

---

Published in Transactions on Machine Learning Research (11/2025)
C.3 Additional Experiments with Non- Non-homogeneous Boundary Conditions
In this section, we add lest case for non-homogeneous boundary conditions using non-homogeneeuu boundary variant of the heat equation to test our proposed framework. Table 20 reports the final relative & errors (L2RE).
PDE Description. boundary conditions:
We use the following non-homogeneoous heat task with mixed Shrichlet/Robin-stte
Unr (xx,1) = (x, () [0, x [0.1]. Ut
Initial condition
u(x, 0) sin(2 + 0.5
Boundary conditions
m(0, 1,
(0(1,))+uell 9(),, h= 4.
Forcing term
f(z, t) = 4m² sin (inx) cos(nt) - (in.2r sin(tt
3
3
5 (6.3x) (6.36))) (3rt) + 3m (6n(6.3 cos(3Tt).
5
Robin data
189T
g(t) 2n cos
cos(Tt) (T))
cos(6.37)
[ (1..6.3) + 1.5] 0.5
50
2n cos(Tt) + 3. cos(6.3T) sin (3rt) + 2.4 2..sin(0.37)sinn33r) + 6.5
Analytical solution
u(x, t) = sin(( cos(It) + sin(6.3Txz) sintt + 0.5 +1
The loss weights, number of features, network widths, sampling strategy, and all other hyperparameters follow the default experimental setup used throughout the paper (See Appendix A). Each result is averaged over five initializations with Optimization Schedule (2).
Training protocol.
Table 20: Non-homogeneous heat. best L2RE results in bold.
Method
L2RE
PINN
2.92e-3
FLS-PINN
3.08e-3
W-PINN
7.01e
RFF
5.12e-3 3
RBF
9.28e 4
RBF-P
3.06e 4
RBA-PINN
3.220
SAFE-NET
2.11e 4
C.4 Spectral Density Computation
In this section, we provide a detailed analysis of the spectral properties of the Hessian matrix for SAFE-NET compared to bascline PINN architectures. We empirically observe that SAFE-NET improves the spectral density in every tested PDE compared to the baselines.
43

---

Published in Transactions on Machine Learning Research (11/2025)
C.4.1 Spectral Density Computation Methodology
The spectral density p()) of a symmetric matrix E _ is defined as:
((A) = Σ 11)
(44)
where A are the eigenvalues of H and 6-) is the Dirac delta function. For neural network optimization, the matrix of interest is the Hessian (0) V2L of the loss function with respect to the network parameters 0.
Computing the full eigendecompposition of the Hessian is computationally prohibitive for neural networks with thousands of parameters. Following the methodology established in Yao et al. (2020) and Golub & Meurant (2009), we employ Stochastic Lanczos Quadrature (SLQ) to efficiently approximate the spectral density. SLQ requires only Hessian-vector products, which can be computed efficiently using automatic differentiation.
The SLQ algorithm approximates the spectral density by:
1 Bp nye
p(A)
of 5 H)U
(45)
where by I) are random Gaussian vectors, and ny is the number of stochastic samples. The quadrature is performed using the Lancsos algorithm, which builds an orthogonal basis for the Krylov subspace suan(u, Hv, H......) and computes the eigenvalues of the resulting tridiagonal matrix.
For our implementation, we utilize the PyHessian library Yao et al. (220)). which provides efficient GPU- accelerated computation of spectral densities for neural networks. We use n, 100 stochastic samples and a Lanezos iteration count of 200 to ensure accurate spectral density cstimation across all experiments.
C.4.2 Experimental Setup for Spectral Analysis
Optimizer Choice: Crucially, we conduct all spectral density experiments using Adam optimizer only, rather than the hybrid Adam + L-BFGS approach used in our performance comparisons. This design choice is essential because L-BFGS inherently performs Hessian preconditioning, which would confound our analysis of the architectural effects on problem conditioning. By using Adam, can isolate the impact of SAFE-NET's feature engineering the fundamental spectral properties of the optimization landscape.
Training Protocol: We train all networks for 100,000 iterations using Adam with initial learning rate = 0.001 and exponential decay factor 0.9 every 2,000 iterations. This extended training schedule ensures that we capture both carly-stage conditioning properties (at 3,000 itcrations) and late-stage spectral characteristics (at 100,000 iterations). All other hyperparameters remain ideutical to the main experimental setup described in Appendix A
Architecture Consistency: To ensure fair comparison, SAFE-NET and all baseline methods maintain identical parameter counts where possible. SAFF-NET uses 128 features as in the main experiments, while baseline methods use their standard configurations, Network weights are initialized using Xavier initialization Glorot & Bengio 2010) across all methods.
Spectral Density Sampling: We compute spectral densities two critical training phases: (1) Early training (3,000 itcrations) to analyze initialization and carly conditioning properties, and (2) End of training (100,000 iterations) to examine how the spectral properties evolve throughout optimization.
C.4.3 Results and Analysis
Figure 10 demonstrates the spectral density comparison between SAFE-NET and standard PINN for the wave equation, showing dramatic improvements in conditioning. The comprehensive results across all PDEs provided visualization examples.
are presented in Tables 21 and 22
Also, Figures 22
26
44

---

Published in Transactions on Machine Learning Research (11/2025)
Table 21: Spectral Properties Comparison: Maximum Eigenvalue and Conditioning Improvements Across Methods and PDEs
Wave Amna
Reaction
Convection
Heat Amax
Burgers
Diffusion
Allen-Cahn
Method
Training Phase
Amaz
Amax
Amax
Amaox
Amna
Early (3k iter)
1.1e+4
4.8e++3
7.2e+3
2.le+3
2.8e+3
1.2e+3
2.5e+3
PINN
End (100k iter)
5.8e+-6
7.le+4
3.9e+-5
8.9e+44
1.6e+44
1.4e+4
2.2e++6
Early (3k iter)
1.8c+44
4.2c-+3
6.5e++3
1.9c+33
2.6c+3
1.c+33
2.3e+3
FLS-PINN
End (100k iter)
4.2c+-6
6.8c+14
2.9c++5
7.8e++4
1.5c+44
1.3c+4
1.9e+66
Early (3k iter)
1.6e+4
3.9e++3
6.1e+3
1.8e+33
2.4e+3
1.0e+33
2.1e+3
W-PINN
End (100k iter)
3.8e+66
6.2e+44
2.6e+-5
7.2e++4
1.4e+-1
1.2e+4
1.7e46
Early (3k iter)
8.9e++3
2. 2e1e+
3.2e+3
9.8+++
1.8e+33
6.5e+2
1.4++3
RBA-PINN
End (100k iter)
1.8e+-5
2.9e++4
1.le+5
3.2e+44
8.9e++3
6.9e+33
2.le+5
Early (3k iter)
1.2e+44
2. 2.8e++3
4.1e+33
1.4e433
2.2e+-3
8.2e+2
1.7e+3
RFF
End (100k iter)
7.8e+-5
4.2e+44
1.8e+55
4.8e4-4
1.1e44
9.8e++3
5.2e+5
Early (3k iter)
1.1ct4
2.6+++3
3.8e+3
1.3c+3
2.00+-3
7.9c+22
1.6e+3
RBF
End (100k iter)
2.9c+-5
3.8c+14
1.5e+5
4.20+44
1.0c+44
9.2c++3
4.1c+5
Early (3k iter)
9.3e+
2.2c++3
3.4e+3
1.1e33
1.7e+
7.0c+
1.5e+3
RIB--P
End (100k iter)
1.5+44
2.8c+44
9.8e
3.1e+44
8.9e+-3
7.8c++3
2.8e++
Early (3k iter) End (100k iter)
1.4e+00
2.1e+
4.5e+
1.3e+
7.8e+11
1.5e+-0
1.le+2
SAFE-NET
4.8e+22
3.5e++3
5.9e+
8.2e++2
8.1e+11
8.9e++1
1.5e+-2
Table 22: Eigenvalue Reduction Factors: SAFE-NET vs. Bascline Methods (End of Training)
vs. SAFE-NET
Wave (x10*)
Reaction (x10')
Convection (x10²)
Heat (x10')
Burgers (x10*)
Diffusion (x10º)
Allen-Cahn
Reduction Factor
(x10*)
PINN
1,2x104
200
6.6x102
1.1x102
2.0x102
1.6x102
1.5x10
FLS-PINN
8.8x
1.9x10'
4.9x102
9.5x
1.9x102
1.5x102
1.3x
W-PINN
7.9x 1
1
104x
1.8x10
1.7x102
1-3
1.1
RBA-PINN
3.8x102
8.3x10
1.9x102
3.9x10
1.1x102
7.8x10
1.4x10"
RFF
16x
1.2x10
3.1x102
5.
1,41022
1.1
3.5x
RBF
6.0x102
1.1x10
2.5x10
5,1x10
1,2x102
100100
2.7x10²
RBF-P
3.1x10
8.0x10
1.71102
3.8x
1.11002
8.8x10
1.9x10
Early Training Analysis: Figures, 10 10 22 a). 23 a) 24 a), 25 a), and 26 a) display spectral density
plots at early training stages (3,000 iterations). Even al this early phase, SAFE-NET demonstrates improved conditioning compared to standard PINNs. The eigenvalue distribution is more concentrated around moderate values, with fewer extremely large cigenvalues that typically cause ill-conditioning. This suggests that SAFE- NET's feature engineering provides inherently better conditioning from initialization, facilitating more stable gradient-basel optimization from the onset of training.
End of Training Analysis:
The spectral density plots at the end of training (Figures 10 10 22 b). b),
24 b). (5(b), and b) reveal even more dramatic conditioning improvements. Specifically for SAFE-NET
and PINN:
Wave and Convection: The largest eigenvalues are reduced by approximately 10 and 10 for the wave and convection lasks respectively, indicating substantial improvement in the condition number of the Hessian matrix.
Heat and Burgers: The largest eigenvalues show reductions of approxinately 10?, demonstrating consistent conditioning benefits across different PDE types.
45

---

Published in Transactions on Machine Learning Research (11/2025)
Eigenvalue Density: Across all problems, SAFE-NET exhibits significant reduction in both the number and density of large eigenvalues in the plots, leading to potentially more favorable optimization landscapes.
REF PRN
PAN
3.0
20'50
33-11
Opanmitie Ac
Ogemmatu
(a) Beginning of Training
(b) End Training
Figure 22: Spectral density plots the beginning and end of training for the wave PDE with RFF and RBF. Compare with Figure I0 demonstrating spectral density for PINN and SAFE-NET
Dexction
Dhaetton
SATE
PAN
SATE
19-
11"04
18
miger.o
(a) Beginning of Training
(b) End of Training
Figure 23: Spectral density plots at the beginning and end of training for the reaction PDE
Contection
Connection
Convection
Connection
100
12:
13'
12/0
3.7
08 02 10
(a) Beginning of Training
(b) End of Training
Figure 24: Spectral density plots at the beginning and end of training for the convection PDE
46

---

Published in Transactions on Machine Learning Research (11/2025)
Polr
Maz:
Hou
Roa
123
10-5
337
(a) Beginning of Training
(b) End of Training
Figure 25: Spectral density plots at the beginning and end of training for the heat PDE
Bargers
Burgers PRN
Burgers
Durgers
142
14-1
307
-07
Beginning of Training
(b) End of Training
Figure 26: Spectral density plots at the beginning and end of training for the burgers PDE
Implications for Optimization: The improved spectral properties could explain SAFE-NET's empirical results despite its simple structure. Well-conditioned Hessian matrices enable stable gradients, faster convergence and reduced sensitivity to learning rate choices and other hyperparameters.
D A Comparison with Large-scale Architectures on Shared Benchmarks
We include the following reported numbers verbatim from the original papers of FNO Li et al. (2020b), PirateNets Wang et al. (2024). PINNsFormer Zhao et al. (2023) as examples large complex models only for the sake of completeness and give readers an at-a-glance sense of scale (parameters, memory, time/cpoch) and accuracy on overlapping PDE families, however, provide direct ranking between them. Note that FNO is an operator -eernin method that learns mappings between function spaces (e.g,, initial condition + solution), whereas SAFE S-FE-NET (and the PINN-style baselines considered in our main text) solve a single PDE instance with specified IC/BC. Formulations, metrics, and datasets therefore differ. The following is intended only to document. the resource scale and the published accuracy on broadly overlapping PDE families.
each method utilizes distinct settings, we do not
FNO Reports results for 1D Burgers and 2D Navier-Stokes (space-time operator learning). Hardware noted by the authors: single NVIDIA V100 GB.
Table 23: FNO on 1D Burgers (relative l2 error al different spatial resolutions s).
Method
8525
512
1024
2048
4096
8192
FNO
0.0149
0.0158
0.0160
0.0146
0.0142
0.0139
Notes. Table reproduced from the paper; parameters, GPU memory, and time/epoch were not reported the Burgers experiment. See Table 24 for Stvier-ttokes resource numbers as reported by the authors.
47

---

Published in Transactions on Machine Learning Research (11/2025)
Table 24: FNO on 2D Navier Stokes (relative & error over different viscosities v and dataset sizes N;
per-epoch time reported by the authors).
Method
Params
Time/epoch
1=10 7, T=500, 1 P=10 T=30, N=1000
y=10 + T=30, N=100000 1=10 TT=20, N=100X
FNO-3D
6,358,538
38.90s
0.0086
0.1918
0.0820
0.1893
FNO-2D
414,517
0.0128
0.1559
0.0973
0.1556
U-Net
29,950.41
48.67
0.0245
0.2051
0.1190
0.1982
Net
7,451,724
47.21
0.0225
0.2253
0.1168
0.2268
ResNet
266,641
78.47s
0.0701
0.2871
0.2311
0.2753
Notes. Reported at 64x64 spatial resolution; FNO-3D convolves in space-time while FNO-2D uses 2D convolations with
RNN in time.
PirateNets has PINN backbone with physics-infornne residual adaptive blocks. The paper emphasizes accuracy comparisons and ablations; it does not tabulate parameter counts, GPU memory, or wall-clock per epoch. Below we list the state-of-the-art test errors the authors report.
Table 25: PirateNets: reported relative le test errors across PDEs (paper's Table 1).
Benchmark
Error (PirateNct)
Params GPU Mem
Time/epoch
Allen-Cahn (1D)
2.24x10-5
Korteweg De Vries (1D)
4.27x10-1
Grey Scott (2D)
3.61x10-=
Ginzburg-Landau (2D)
1.49x 10
Lid-driven cavity (2D)
4.21x10-
Notos. Architecture details (o.g., width) and training pipelines provided, but metrics not tabulated.
PINNsFormer is transformer-style PINN variant. The authors report parameter counts and training overhead (V100), and test errors overlapping 1D PDEs.
Table 26: PINNsFormer: model size and training overhead (Appendix Table 4-5 in the paper).
Model
Params
GPU Mem (MiB)
Time/epoch (s)
PINNsFormer (pseudo-seq. length k=5)
454,000
2,827
2.34
Notes. Reported single NVIDIA Tesla V100; overheads shown for k=5.
Table 27: PINNsFormer: reported lest errors on 1D PDEs used widely PINN literature.
PDE (dimension)
Metric (paper)
Error
Params
Time/epoch / GPU Mem
Convection (1D)
rRMSE (ss rel. f2)
0.027
454k
2.34s 2,827 MiB
Reaction (ID)
rRMSE (st rel. 62)
0.030
454k
2.34s / 2,827 MiB
Notes. Errors taken directly from the paper's main results tables; rRMSE the paper's standard relative metric. The paction/convectiin formulations and sampling follow setups specified Zhao all 2023
48

---

Published in Transactions on Machine Learning Research (11/2025)
E A Discussion on Relevant Ideas in Different Structures
In this section we discuss three different frameworks that leverage Fourier structure for PDEs in different ways; SAFE-NET, FNO Li et al. (2020b), and EPGP Härkönen et al. (2023). Although all three exploit Fourier structure, they do so for different learning objectives and with complementary inductive biases.
Problem formulation and learning objective. SAFE-NET a physics-inlormee solver for a single PDE instance: we fit ug by minimizing a residual-basee PINN loss with collocation points for the PDE, BCs and (if present) ICs; no paired input-output simulation dataset required. FNO Li et al. (2020b), by contrast, learns ansolution-invarian solution operator G + from many paired fields I(a 2 Ti once trained it solves new instances in one forward pass and even on unseen grids (zero-shot super-resolution, EPGP from Härkönen et al. (023) is probabilistic: it places a GP prior whose sample paths satisfy a linear, constant-cocfficien PDE exactly; conditioning on scattered observations (valucs and/or derivatives) yields a posterior that remains within the PDE solution set.
Spectral parametcrization. All three approaches loverage Fourier structure, but in different ways.
SAFE-NET uses a finite, trainable tensor-product Fourier feature map in space (and time), with sine and cosine cross terms and learned amplitudes. Both the harmonically-initaann and amplitudes are optimized jointly with the network via the physics loss. This couples spectral adaptivity to the governing equation and moundayboundary condition information embedded in the loss.
frequencies
FNO parameterizes convolutions in Fourier space: the FFT of the feature field multiplied by learnable complex matrices on a fixed set of retained modes and transformed back by an inverse FFT. In other words, each Fourier layer applies (x) -F- (R)) + We((z), where R is a learnable complex tensor restricted to low modes and W local lincar map; frequencies themselves remain lixed and are not learned; instead, (problen-dependent)) truncation kmax selects the modes that are updated, while higher modes are filtered out Li et al. (2020b). This yields global, resolution-invarian layers well suited to operator learning.
EPGP in Härkönen et al. (0223) expressos solutions through the Ehrenpreis Pallaooo representa- lion (a spectral integral over the PDE's characteristic variety) as the following kernel
kEPGP (0,2') = -1, JV-IR= (x 2) V (z', 2) -1=11//22
(46)
with ((x
= = = Dj(2, 2) exp(x, 2). A Monte-Carlo variant (S-EPGP) replaces the integral by a finito sum of trainable planc waves, so the relevant spectral locations are learned within the GP kernel.
Data/physics requirements and scope. SAFE-NET relies on the PDE residual (and BC/IC) and targets cow.toemodderaae-dimeniioo bases provide a good inductive bias. FNO excels when sufficiently many input output pairs of fields are available; its strength is generalization across instances and grids rather than solving a single instance from physics alone. EPGP is model-basee and probabilistic: it assumes linear PDEs with constant coefficients and incorporates data at arbitrary locations (values and/or derivatives).
PDEs (including smooth and mildly non-smooth regimes), where Fourier
Boundary/initial conditions and domain assumptions. SAFE-NET feeds BC/IC conditions, if available, directly as domain features (c.g., sinc scries that satisfy homogencous Dirichlet BCs) and trains end-to-end with the PINN loss. FNO typically assumes periodic structure within each Fourier layer (the convolutional kernel is defined for periodic functions), although practical variants can handle nonperiodic domains with padding or localized corrections. EPGP does not require prescribing BCs BCs/ICss define the prior; when BC/IC observations are available, they are simply assimilated as GP conditioning data, and the posterior remains in the PDE solution set.
49

---

Published in Transactions on Machine Learning Research (11/2025)
Computational aspects. Each of these architectures solve different problems they should not be compared directly based on computational cost. However, for the sake of completion, provide computational footprint details for each method. SAFE-NET uses a single hidden layer and a small trainable feature module, giving low parameter counts and fast epochs suitable for residual-base bassed training. FNO stacks multiple FFT-based layers; training scales as O(Nn log n) but inference on new grids is extremely fast. EPGP shifts cost to kernel construction O(n²) (or Onm²² for sparse variant) while providing calibrated uncertainty. With sparse/Monte Carlo kernels, the dominant cost is one-time covariance construction.
Complementarity. FNO works well when abundant multi-instance data are available, offering orders- of-magnitude speed-ups at deployment and results generalizable within the parametrized family of PDEs. SAFE-NET is best when the PDE (best for PDEs with moderately smooth behavior) is known but training data are searce, computational budget limited, or when we seek both accuracy and efficiency together. EPGP delivers principled uncertainty for linear PDEs and fuses scattered measurements seamlessly.
50


Reviews:
[{"review_id": "OYwPZuDVeK", "reviewer": "Reviewer_9qG6", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: The proposed model is very simple and clean, yet it achieved great numerical performance over many PDEs. The author also did a comprehensive comparison study with different feature engineering approaches over several benchmark PDEs.\n\nWeakness:\nA major problem of this approach is its scalability in dimension. Even if for a 1D time-dependent PDE, the number of basis N in (9) can easily grow beyond 100 with only ten quadrature points per dimension in the 2D Fourier transformer. \n\nThe preconditioning analysis in Section 6 is built up unrealistic assumptions, making its justification of the better conditioning weak. Please see details in requested changes.\n\nThere are quite many typos and inconsistent notations throughout the paper, including but not limited to the ones listed below. The authors should make a further screening to avoid them.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposed a single layer network coupled adaptive fourier feature maps to improve the condition of PINN optimization. It is shown that such small network can attain decent accuracy and speed for most linear PDE problems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Some of the testing cases are oversimplified as most PDEs are equipped with initial conditions of only one frequency and homogeneous boundary condition. Appendix A.1.2. Most equations are tested with very simple trigonometric initial conditions and homogeneous boundary conditions. Such oversimplification can not fully demonstrate the effectiveness of the proposed SAFE-NET. \n\nThe preconditioning analysis in section 6 may not be suitable as it used NTK as the main tool, where it is assumed that the number of model parameters p goes to infinity. This assumption may not be consistent with small size of SAFE-NET, which is one of the main advatange of SAFE-NET over other feature engineering approaches. Furthermore, the preconditioning analysis did not consider adaptive frequencies w and \\lambda as model parameters in (11), but directly assuming that \\phi^{(i)}_j are orthonormal. This assumption is questionable as the learned frequcies might not be integer multiples of \\pi.\n\n\nTypos:\nbelow eqn: Where -> where\nabove (13): \\Theta_n -> \\Theta_u\nMisspelled of Convection in Table 3.\n\nNotations:\nThe formulation (1) is confusing for time-dependent problems; For example, the 2nd equation should be defined on \\partial\\Omega \\times (temporal domain).\nNotations in (3) and (4) are not consistent.\nPage 5: Frequency notations are not consistent in Figure 2 and the discussion below it.\nSuperscript of frequncies are somehow missing in (9).", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "YWosjUu25e", "reviewer": "Reviewer_Pn2M", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The strength of the paper is the ablation study done on the methods of training as can be seen in Appendices C and A.3. \n\nThe key weaknesses of the paper are the misplaced choices of demonstration and lack of rigour in making the point about ill-conditioning.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The key idea of the work is train PINNs using simple model which hand-crafts a functional transform of the inputs and then passes it on to a single layer of activation. The hand-crafted features are a linear combination of Fourier-like modes in space and temporal variables with learnable frequencies. \n\nThe experimental evidence is very thin in terms of covering conventional hard cases although for the simple exemplars chosen the authors have given extensive ablation studies for the algorithmic and the architectural choices.\n\nAnd even among the chosen simple tests, the given data is ambiguous --- for 50\\% of the PDEs i.e Diffusion, Allen-Cahn, Burgers and Navier-Stokes the data in Table-3 says that SAFE-NET is not the one with the lowest error. So it seems that the author's own data undermines their proposed model! \n\nBut then again these could be just strange artefacts of the unexplained capping of the models to just one layer of activation - this choice is possibly fundamentally affecting the results but that has never been highlighted or explained!\n\nThere is some attempt to justify the idea in terms of saying that \"SAFE-NET enjoys an implicit preconditioning effect that is global\" - but there is nothing rigorous about this idea visible in the draft - not even a formally defined conjecture. \n\nAs it stands, this paper needs a significant and almost total overhaul to be accepted.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The paper needs to be entirely rewritten/recalibrated to include non-trivial PDE examples. The demonstration choices made (Table 2) are not at all well-chosen to serve as an indicator of performance. Of the 8 PDEs chosen, none have 3 or more variables and 6 are having only one or 2 variables and these are not even known to be classically hard PDEs. \n\n  So, firstly we need to see evidence of performance in high number of variables. PDEs like Burger, anisotropic linear Fokker-Planck, HJB and heat are some non-trivial examples where exactly solvable cases are available in every dimension. Hence these can be good test cases to benchmark against --- particularly the non-linear ones.  \n\n- One can see the kind of examples that were chosen for demonstration in this work, https://arxiv.org/abs/2407.20741. The examples here not only consist of 4 variable PDEs but also consider classically hard PDEs for even 2 variable PDEs like the 3-solitonic solution for KdV and the methods therein were shown to supersede the erstwhile best known solution for that in terms of needing smaller models. Such a demonstration is needed to make SAFE-NET a credible competitor and for that these same benchmarks can be tested against. \n\n- Also for the simple one and 2 variable PDEs chosen in this work, the boundary and initial conditions are simple enough that they could be incorporated into the model as shown in the demonstrations in the above paper. So, with SAFE-NET, we definitely need this study whether SAFE-NET with inclusion of boundary and initial conditions in the model supersedes a standard PINN with the same inclusions in the model.\n\n- The point about the ill-conditioning needs to be much more rigorously said. If the authors want to claim that \"SAFE-NET enjoys an implicit preconditioning effect that is global\" then this needs to be said as a mathematical conjecture - currently even the notion of preconditioning as in this context has not been defined. In Appendix C.3.1 there is a discussion on the spectral density, but its not clear as to what is the formal conjecture about the spectral density of the Hessian of SAFE-NET that they want to make. There is a need for a lot more clear conjectures and showing how the evidence aligns with it. \n\n- Lastly, its entirely unclear why the authors are stuck to using only one layer of activation. Maybe SAFE-NET performs great at high-depth and its not clear why these obvious experiments were not done. Or, if the authors want to claim that among all depth-2 models for PDE solving, SAFE-NET is the optimal choice, then that needs to be said that way and the whole paper needs to be positioned that way. Currently this capping of the depth to 2 looks entirely artificial. \n\n\n\nMinor Issues,\n\n- the definition of the tangent kernel in Definition 6.1 seems to have a typo wih its evaluation points.\n\n- Section 5.1 in its opening line attempts to classify the PINN baselines in 2 categories but that organization is not maintained when the actual review is done in Appendix A.2. This is confusing writing and needs looking into. When defining RBF-PINN, the variable $c_i$ has not been defined properly. \n\n- For the one variable setup, a classic hard example to test would be the energy eigenfunctions of the quantum harmonic oscillator. At quantum numbers $3$ or more, these are strong targets for contest between vanilla PINNs  and any other modification.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "eeS5Bom4bM", "reviewer": "Reviewer_LL46", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n\n1. The experiments are very fully described.\n2. Making many ablations w.r.t.optimization schedules and an NTK-based theoretical discussion about loss landscapes, which also give additional insights into problems with PINN-training\n\nWeaknesses:\n\n1. While some justification is given as to why no large scale experiments are being performed, I would be intereed in at least 1-5 training runs.\n2. While some intuition is given about the features, I cannot say that the deep connection between choosing the features precisely this way and the theory of differential equations.\n\nQuestions\n\n1. I am confused about the notation in your first paragraph. I a potentially existing time part of the domain $\\\\Omega$? If so, the initial condition should probably just be evaluated at $\\\\partial\\\\Omega$. Can you clarify?\n2. You are considering PINNS \"just\" for solving PDEs. They are also a regression tool in the context of PDE, in case data is given. Perhaps you might want to clarify your stance on this distinction. Does your approach also work in this regression setting? In particular, you seem to assume knowledge about the Fourier transform of initial and boundary conditions. This might or might not be a reasonable assumption, depending on circumstances. Can you comment here?\n3.  You are learning a (non-regular) grid of Fourier frequencies. How does this scale to high input dimensions? Did you try experiments for higher dimensions? (Usually, an advantage of PINNs over methods from scientific computing is, that PINNs are less vulnerable against the curse of dimensionality. Here, it seems, you get a strong curse of dimensionality.) The question holds particularly since your experiments use at most 2 input dimensions (Table 2, wven though if I understand correctly, having \"time dependent\" here means you use an additional input dimension?). Does your approach scale to 3 or 4 dimensions? How long would computations take? (Yes, I read your comment at the end of Section 3, but can you say more here?)\n4. Fourier based theory for differential equations is most relevant for linear differential equations. Can you comment on differences in the performance of your approach between linear, close-to-linear, and highly non-linear PDEs? See also Figure 3, where Safe-Net performes particularly well for the linear equations and at best mediocre for non-linear equations.\n5. Are the learned frequencies in any way related to the frequencies in the characteristic varieties of the linear PDEs or to the Fourier frequencies of the initial conditions?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper constructs specific features suitable for differential equations, such that PINNs are significantly easier to train using these features. These feature use Fourier transforms and Fourier series, amongst others. Careful experiments are conducted on small neural networks, which show that the method SAFE-NET, introduced in this paper, is superior in quality to baseline models, in particular for linear equations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. You might want to check your placement of figures. Sometimes they appear on seemingly random positions of the page.\n2. The text in the images used in the figures is much too small to be readable. Please enlarge it at least to the size of footnotes.\n3. Introduction: \"PDE solutions have well-understood mathematical structures rooted in Fourier analysis, differential geometry, and physics principles\" I think this sentence might profit from references into the literature, or making an entire paragraph out of it.\n4. On page 6, the systematic feature extraction method is described really vague. How does this work?\n5. Page 6 comments on \"prioritizing\" certain frequencies for homogeneous Dirichlet boundary conditions. How does this work in detail? Can you also comment on the connection of this to Fourier series?\n6. While I get that comparisons against mall models is reasonale here, can you comment on accuracy, traning times and memory for larger models, e.g. by citing their papers?\n7. In the first, unnumbered, equation on page 15, you have a limit $p\\\\to\\\\infty$, but no $p$ appears afterwards in the equation (ok, it is the number of parameters, but this is not obvious locally.).\n8. In appendix A.1.4, using i as a running variable is somewhat confusing, since $i=\\\\sqrt{-1}$ also makes sense.\n9. In appenxis B.3, there is a typi $x in ...$ instad of $x\\\\in ...$.\n10. Some relations hold to various neural operator approaches, e.g. FNO, which you are citing. Can you comment more on similarities and differences?\n11. A similar approach called EPGP also trained Fourier frequencies and amplitudes for Gaussian processes, but not in a grid. Can you comment on similarities/differences?\n12. The paper \"Sifan Wang, Xinling Yu, and Paris Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective. Journal of Computational Physics\" appears twice in the bibliography.\n13. The paper \"Attention is all you need\" is not from 2023, and it is not just a preprint.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)