Operator-informed initialization for Fourier features physics-informed neural networks

Juan Molina; Paris Perdikaris; Mircea Petrache; Matías Courdurier; Francisco Sahli Costabal

cs.LGOct 2, 2026

Operator-informed initialization for Fourier features physics-informed neural networks

Authors: Juan Molina, Paris Perdikaris, Mircea Petrache, Matías Courdurier, Francisco Sahli Costabal

Organizations: National Center for Artificial Intelligence, Santiago, Chile · Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, PA, USA · Department of Mathematics, Pontificia Universidad Católica de Chile, Santiago, Chile · Institute for Mathematical and Computational Engineering, Pontificia Universidad Católica de Chile, Santiago, Chile ·

Abstract

Physics-Informed Neural Networks (PINNs) typically exhibit spectral bias, where some frequencies of the target function converge more slowly than others. In this work, we analyze the training dynamics of Fourier Feature PINNs in the Neural Tangent Kernel regime to address this limitation. We derive an explicit evolution equation to estimate the residual error in the frequency domain, demonstrating that the convergence rate of specific frequencies is primarily governed by the product of the differential operator's symbol and the spectral density of the initialization weights. Leveraging this theoretical insight, we propose an informative initialization strategy that tailors the initial weight distribution to the specific PDE being solved. With this method, we can diminish the operator-induced spectral bias, balancing the convergence rates across the frequency spectrum and achieving better prediction accuracy. Numerical experiments on linear and nonlinear partial differential equations confirm that this initialization strategy improves learning dynamics and approximation accuracy across frequencies compared to standard initialization methods, with no additional training cost.

Figures & tables

Panel of Figure 1 — Figure 1: Pipeline of the proposed PDE-informed initialization. The governing PDE determines the principal symbol $p_{\mathrm{pde}}$ , which is used to sample the Fourier frequencies $W$ through the proposed initialization density. The sampled frequencies define the Fourier Feature embedding of the PINN, whose parameters are trained by minimizing the PDE, boundary, and initial-condition residual losses.

Figure 2: Frequency-domain (modulus) error dynamics for three PDEs solved with PINNs. From top to bottom: the 2D Helmholtz equation

u_{xx}+u_{yy}+k^{2}u=0

, with

k=10

, and the 1D transport equation

u_{\tau}+cu_{x}=0

with

c=1

and

c=80

. The corresponding polynomials are

p_{\mathrm{pde}}=(4\pi^{2}(\xi_{x}^{2}+\xi_{y}^{2})-k^{2})^{2}

for the Helmholtz equation and

p_{\mathrm{pde}}=4\pi^{2}(\xi_{\tau}+c\xi_{x})^{2}

for the transport equation. The results show that the largest residual errors tend to persist in frequency regions where

p_{\mathrm{pde}}

is small or vanishes, consistently with the weak interior control predicted by equation ( 27 ).

Figure 3: Visualization of the exact solutions listed in Table 1 . The examples cover a diverse range of spatial structures and spectral characteristics representative of the benchmark PDEs considered in this work.

PDE	$p_{\mathrm{pde}}$	Exact solutions
Transport $(c=1)$ $u_{\tau}+c\,u_{x}=0$	$(\xi_{\tau}+c\xi_{x})^{2}$	i) Multiscale oscillatory: $u(\tau,x)=\cos(11\pi(x-c\tau))+\cos(\pi(x-c\tau))$ in $\Omega_{3}$ ii) Localized Gaussian: $u(\tau,x)=e^{-100(x-c\tau-0.2)^{2}}$ in $\Omega_{2}$
Heat $(\alpha=0.1)$ $u_{\tau}-\alpha u_{xx}=0$	$\xi_{\tau}^{2}+4\pi^{2}\alpha^{2}\xi_{x}^{4}$	i) Single-mode: $u(\tau,x)=e^{-\pi^{2}\alpha\tau}\sin(\pi x)$ in $\Omega_{3}$ ii) Multiscale: $u(\tau,x)=e^{-\pi^{2}\alpha\tau}\sin(\pi x)+\frac{1}{11}e^{-121\pi^{2}\alpha\tau}\sin(11\pi x)$ in $\Omega_{2}$
Wave $(c=1)$ $u_{\tau\tau}-c^{2}u_{xx}=0$	$(\xi_{\tau}^{2}-c^{2}\xi_{x}^{2})^{2}$	i) Multiscale oscillatory: $u(\tau,x)=\cos(\pi\tau)\cos(\pi x)+\cos(11\pi\tau)\cos(11\pi x)$ in $\Omega_{2}$ ii) Traveling Gaussian: $u(\tau,x)=\frac{1}{2}\left(e^{-50(x-c\tau-\frac{1}{2})^{2}}+e^{-50(x+c\tau-\frac{1}{2})^{2}}\right)$ in $\Omega_{3}$
Laplace $-\Delta u=0$	$(\xi_{x}^{2}+\xi_{y}^{2})^{2}$	i) Oscillatory-hyperbolic: $u(x,y)=\sin(\pi x)\sinh(\pi y)$ in $\Omega_{1}$
Helmholtz ( $k=10$ ) $\Delta u+k^{2}u=0$	$(4\pi^{2}(\xi_{x}^{2}+\xi_{y}^{2})-k^{2})^{2}$	i) Multiscale oscillatory: $u(x,y)=\cos(kx)+\cos(ky)+\cos\left(\frac{\sqrt{2}k}{2}(x+y)\right)$ in $\Omega_{1}$
Nonlinear Laplace $\Delta u+u^{3}=0$	$(\xi_{x}^{2}+\xi_{y}^{2})^{2}$	$u(x,y)=\dfrac{1}{1+x^{2}+y^{2}}$ in $\Omega_{2}$

Table 1: Benchmark PDEs used for validation. The exact solutions are selected to exhibit diverse spectral characteristics representative of each PDE. The computational domains are

\Omega_{1}=[-1,1]^{2}

\Omega_{2}=[0,1]^{2}

, and

\Omega_{3}=[0,1]\times[-1,1]

. For simplicity, multiplicative constants in the polynomial

p_{\mathrm{pde}}

are omitted.

example	$\eta_{\mathrm{GI}}$	$\sigma_{\boldsymbol{w}}^{\mathrm{GI}}$	GI range	$\eta_{\mathrm{PI}}$	$\sigma_{\boldsymbol{w}}^{\mathrm{PI}}$	PI range
Helmholtz Multiscale Osc.	$5\cdot 10^{-7}$	$1.42$	$[0.3,3]$	$5\cdot 10^{-4}$	$0.028$	$[0.005,0.5]$
Laplace Osc.–Hyperbolic	$8\cdot 10^{-7}$	$1.4$	$[0.5,5]$	$2\cdot 10^{-6}$	$7.3\cdot 10^{-5}$	$[5\cdot 10^{-5},5\cdot 10^{-4}]$
Trans. Local. Gaussian	$10^{-5}$	$3.0$	$[0.1,5]$	$2\cdot 10^{-3}$	$0.62$	$[5\cdot 10^{-2},5]$
Trans. Multiscale Osc.	$10^{-6}$	$2.1$	$[1,3]$	$8\cdot 10^{-5}$	$2.6$	$[0.1,10]$
Heat Single-Mode	$5\cdot 10^{-5}$	$0.34$	$[0.1,5]$	$5\cdot 10^{-5}$	$7.7\cdot 10^{-3}$	$[10^{-4},10^{-2}]$
Heat Multiscale	$6\cdot 10^{-5}$	$0.25$	$[0.1,1]$	$6\cdot 10^{-4}$	$0.78$	$[0.1,1]$

Table 2: Learning rates used in the experiments and Optuna search intervals for the proposed PDE-informed initialization. The examples correspond to Figure 5 .

Figure 4: Examples illustrating the decoupling between empirical training loss and prediction accuracy. The first row reports the training loss, while the second row shows the relative

L^{2}

prediction error. In the first column, the Gaussian initialization achieves lower training loss but a larger prediction error. In the second column, the PDE-informed initialization achieves lower training loss, whereas the Gaussian initialization yields the more accurate solution.

Figure 5: GD training dynamics for linear and nonlinear PDEs under the standard Gaussian and the proposed PDE-informed initialization. The odd rows report the training loss, while the even rows show the relative $L^{2}$ prediction error. The first two pairs of rows correspond to linear PDEs. The first pair includes the Helmholtz, transport, heat, and wave equations, while the second pair includes the Laplace, transport, heat, and wave equations. For the transport, heat, and wave equations, the two pairs use different exact solutions, thereby assessing the initialization across distinct spectral contents for the same governing PDE. The last pair of rows presents four nonlinear examples: a manufactured nonlinear Laplace equation, the Burgers equation, the Sine–Gordon equation, and the nonlinear Klein–Gordon equation. For the Laplace equation, both initialization strategies exhibit similar dynamics. In contrast, the proposed PDE-informed initialization yields moderate improvements for the heat equation and more pronounced gains for the Helmholtz, transport, and wave equations. For the nonlinear problems, the proposed initialization generally improves prediction accuracy and stability across random seeds.

Figure 6: Evolution of the relative ring-wise spectral error during training for four representative PDEs under the proposed PDE-informed and standard Gaussian initializations. At each iteration, the prediction error is transformed into the Fourier domain, and its relative energy is computed over five concentric radial-frequency bands with logarithmically spaced radii, excluding the zero-frequency mode. The vertical axis orders these bands by increasing radial frequency, from the innermost low-frequency band to the outermost high-frequency band, while the horizontal axis denotes the training iteration. Thus, the vertical coordinate represents an ordered frequency band rather than the actual frequency magnitude or ring radius. For each PDE, the PDE-informed and Gaussian heatmaps share the same color scale, allowing a direct comparison of their spectral learning dynamics. The PDE-informed initialization generally produces a more uniform reduction of the spectral error across frequency bands, with larger differences for problems whose low-symbol regions exhibit stronger anisotropy or nontrivial characteristic geometries.

example	$\eta_{\mathrm{GI}}$	$\sigma_{\boldsymbol{w}}^{\mathrm{GI}}$	GI range	$\eta_{\mathrm{PI}}$	$\sigma_{\boldsymbol{w}}^{\mathrm{PI}}$	PI range
Helmholtz Multiscale Osc.	$5\cdot 10^{-3}$	$0.97$	$[0.3,3]$	$1\cdot 10^{-3}$	$0.71$	$[0.01,1]$
Trans. Local. Gaussian	$5\cdot 10^{-4}$	$2.1$	$[0.1,10]$	$5\cdot 10^{-4}$	$0.096$	$[10^{-4},1]$
Heat Single-Mode	$10^{-4}$	$0.27$	$[0.1,10]$	$10^{-4}$	$0.086$	$[10^{-5},0.1]$
Wave Multiscale Osc.	$10^{-3}$	$2.6$	$[0.5,10]$	$7\cdot 10^{-4}$	$53.5$	$[0.01,100]$
Nonlinear Cubic	$7\cdot 10^{-5}$	$0.095$	$[0.01,10]$	$10^{-4}$	$0.03$	$[0.005,0.5]$
Nonlinear Burgers	$10^{-4}$	$0.11$	$[0.05,0.5]$	$10^{-4}$	$7.78$	$[1,10]$

Table 3: Learning rates, selected initialization scales, and Optuna search intervals used in the experiments. The examples correspond to Figure 7 .

Figure 7: ADAM training dynamics for linear and nonlinear PDEs under the standard Gaussian and the proposed PDE-informed initialization. The odd rows report the training loss, while the even rows show the relative $L^{2}$ prediction error. The linear experiments include the Helmholtz, transport, heat, and wave equations, while the remaining rows present four nonlinear PDEs: a manufactured nonlinear Laplace equation, the Burgers equation, the Sine–Gordon equation, and the nonlinear Klein–Gordon equation. Consistent with the GD experiments, the proposed initialization improves the optimization dynamics and prediction accuracy for the linear problems, with moderate gains for the heat equation and more pronounced improvements for the Helmholtz, transport, and wave equations. Similar improvements are observed for most of the nonlinear problems, with generally faster convergence, lower prediction errors, and more stable behavior across random seeds.

Figure 8: Multilayer experiments under GD and ADAM optimization. The first two rows show the gradient descent training dynamics, while the third and fourth rows correspond to ADAM. The odd rows report the training loss, whereas the even rows show the relative $L^{2}$ prediction error. Both optimizers are evaluated on the transport, heat, nonlinear Laplace (cubic), and nonlinear Klein–Gordon equations. Consistent with the two-layer experiments, the proposed PDE-informed initialization improves both the optimization dynamics and predictive accuracy under both optimization methods. The improvements are particularly pronounced for the transport equation, while more moderate but still noticeable gains are observed for the heat equation. Similar improvements are also observed for the two nonlinear problems, demonstrating that the benefits of the proposed initialization extend to deeper neural network architectures and adaptive optimization methods.

Figure 9: Training dynamics on the JAX-PI benchmarks using the recommended SOTA configuration. The first row shows the training loss, while the second row reports the relative $L^{2}$ prediction error. The columns correspond to the advection, Allen–Cahn, and Kuramoto–Sivashinsky equations, respectively. The blue curves denote the proposed PDE-informed initialization (using ( 34 ) for the linear problem and ( 32 ) with $\lambda=0.5$ for the nonlinear problems), while the orange curves correspond to the standard Gaussian initialization. Solid lines represent the median over ten random seeds and the shaded regions indicate the range between the best and worst runs.

Figure 10: Training dynamics obtained with the proposed PDE-informed initialization under the simplified baseline configuration. The layout is identical to Figure 9 . The results demonstrate that the baseline configuration achieves optimization dynamics comparable to those of the full SOTA pipeline, with the Allen–Cahn equation exhibiting noticeably faster convergence during the early stages of training.

Figure 11: Qualitative comparison of the reconstructed solutions obtained with the SOTA configuration for Gaussian initialization and baseline configuration for PDE-Informed initialization. Each row corresponds to one benchmark equation (advection, Allen–Cahn, and Kuramoto–Sivashinsky). The columns show the ground-truth solution, the prediction obtained with the proposed PDE-informed initialization, the prediction obtained with the standard Gaussian initialization, the absolute error of the PDE-informed initialization, and the absolute error of the Gaussian initialization. The predictions produced by both methods closely match the ground truth, while the proposed initialization generally yields lower reconstruction errors, particularly for the advection equation.

Appendix figures & tables3 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 12: Effect of using an inaccurate PDE symbol in the PDE-informed initialization for the Transport – Localized Gaussian example of Table 1 , whose true transport velocity is

c=1

. The left panel shows the training loss and the right panel the relative

L^{2}

prediction error. The standard Gaussian initialization is shown in orange, while the PDE-informed initialization constructed using the correct velocity

c=1

is shown in blue. PDE-informed initializations constructed using the incorrect values

c=2

c=5

, and

c=10

are shown using progressively lighter shades of green. The PDE used during training always has

c=1

; only the polynomial used to construct the initialization density is modified. Initialization based on the correct PDE symbol yields the best dynamics among the tested alternatives, while performance progressively deteriorates as the assumed transport velocity deviates further from its true value.

Ablation Settings				Performance
Time Period	Grad Norm	RWF	Causal Rel.	L2 Error	Run Time (min)
✓	✓	✓	✓	3.91e-04	22.01
$\times$	✓	✓	$\times$	4.89e-05	22.9
$\times$	✓	$\times$	✓	6.39e-05	21.94
$\times$	$\times$	✓	✓	5.75e-05	21.3
$\times$	✓	✓	✓	5.36e-05	22.02

Appendix

Table 4: Ablation study for the Advection equation.

Ablation Settings			Performance
Grad Norm	RWF	Causal Rel.	L2 Error	Run Time (min)
✓	✓	✓	2.782e-04	43.3
✓	$\times$	$\times$	7.924e-05	42.6
✓	✓	$\times$	1.475e-04	44.6
✓	$\times$	✓	1.153e-04	40.8
$\times$	✓	✓	5.197e-04	40.5

Appendix

Table 5: Ablation study for the Allen Cahn equation.

Physics-Informed ML

Physics-informed machine learning combines data-driven models with physical laws to solve, predict, and infer scientific systems.

See all 833 papers

Feb 16, 2026

Pseudo-differential-enhanced physics-informed neural networks

We present pseudo-differential enhanced physics-informed neural networks (PINNs), an extension of gradient enhancement but in Fourier space. Gradient enhancement of PINNs dictates that the PDE residual is taken to a higher differential order than prescribed by the PDE, added to the objective as an augmented term in order to improve training and overall learning fidelity. We propose the same procedure after application via Fourier transforms, since differentiating in Fourier space is multiplication with the Fourier wavenumber under suitable decay. Our methods are fast and efficient. Our methods oftentimes achieve superior PINN versus numerical error in fewer training iterations, potentially pair well with few/fixed samples in collocation, and can achieve breakthrough rather than gradual effects on the error. Moreover, our methods are suitable for fractional derivatives. We establish that our methods, due to the dynamical effects, improve spectral eigenvalue decay of the neural tangent kernel (NTK), and so our methods contribute towards the learning of high frequencies in early training, mitigating the effects of frequency bias up to the polynomial order and possibly greater with smooth activations. In particular, our primary contribution can be characterized as gradient enhancement affects spectral bias, and we specialize to Fourier space, although our theory does not prove the ambient case exactly due to a result with the Plancherel theorem on the descent paths. Our methods accommodate advanced techniques in PINNs, such as Fourier feature embeddings. A pitfall of discrete Fourier transforms via the Fast Fourier Transform (FFT) is mesh subjugation, and so we demonstrate compatibility of our methods for greater mesh flexibility and invariance on alternative Euclidean and non-Euclidean domains via Monte Carlo methods and otherwise, although possessing trade-offs of their own.

Andrew Gracyk

Department of Mathematics, Purdue University, West Lafayette, IN 47907, United States

PDE Solving Spectral Bias Physics-Informed ML

Ask PDF

Aug 6, 2026

Alternating Levenberg-Marquardt Training of Physics-Informed Neural Networks with Fourier-Enhanced Features

Physics-informed neural networks (PINNs) often fail to accurately resolve partial differential equations (PDEs) with high-frequency or multi-scale solutions, as well as strongly nonlinear problems. Two factors underlie this difficulty: spectral bias, the tendency of neural networks to underfit high-frequency features; and representation-coefficient coupling, the entanglement of representation learning and coefficient fitting within a single nonconvex optimization objective. In this work, we propose the Fourier-enhanced alternating Levenberg--Marquardt PINN (FALM-PINN), an optimization framework that decouples representation learning from coefficient fitting. The upper-level problem learns a Fourier-enhanced basis that enriches the latent space with high-frequency components, while the lower-level problem resolves the coupling by fitting the projection coefficients on this basis, solving a nonlinear least-squares problem with the Levenberg--Marquardt algorithm. The framework applies to general nonlinear and coupled PDE systems, and reduces to a single-step convex optimization problem for linear PDEs. We prove global convergence of the alternating training scheme in both cases. Numerical examples on multiple challenging high-frequency and nonlinear PDEs show that FALM-PINN achieves relative

L^2

errors up to two orders of magnitude lower than state-of-the-art baselines.

Yulun Wu, Matthieu Barreau, Miguel Aguiar +1

Division of Decision and Control Systems, School of Electrical Engineering and Computer Science, KTH Royal Institute of Technology, Stockholm, 100 44, Sweden

Neural PDE Solvers Spectral Bias Physics-Informed ML

Ask PDF

Feb 22, 2026

Spectral bias in physics-informed and operator learning: Analysis and mitigation guidelines

Solving partial differential equations (PDEs) by neural networks as well as Kolmogorov-Arnold Networks (KANs), including physics-informed neural networks (PINNs), physics-informed KANs (PIKANs), and neural operators, are known to exhibit spectral bias, whereby low-frequency components of the solution are learned significantly faster than high-frequency modes. While spectral bias is often treated as an intrinsic representational limitation of neural architectures, its interaction with optimization dynamics and physics-based loss formulations remains poorly understood. In this work, we provide a systematic investigation of spectral bias in physics-informed and operator learning frameworks, with emphasis on the coupled roles of network architecture, activation functions, loss design, and optimization strategy. We quantify spectral bias through frequency-resolved error metrics, Barron-norm diagnostics, and higher-order statistical moments, enabling a unified analysis across elliptic, hyperbolic, and dispersive PDEs. Through diverse benchmark problems, including the Korteweg-de Vries, wave and steady-state diffusion-reaction equations, turbulent flow reconstruction, and earthquake dynamics, we demonstrate that spectral bias is not simply representational but fundamentally dynamical. In particular, second-order optimization methods substantially alter the spectral learning order, enabling earlier and more accurate recovery of high-frequency modes for all PDE types. For neural operators, we further show that spectral bias is dependent on the neural operator architecture and can also be effectively mitigated through spectral-aware loss formulations without increasing the inference cost.

Siavash Khodakarami, Vivek Oommen, Nazanin Ahmadi Daryakenari +2

Division of Applied Mathematics, Brown University, Providence, RI 02912, USA · School of Engineering, Brown University, Providence, RI 02912, USA

Spectral Bias Physics-Informed ML PDE Operator Learning

Ask PDF

Abstract

Physics-Informed ML

Explore similar work

Pseudo-differential-enhanced physics-informed neural networks

Alternating Levenberg-Marquardt Training of Physics-Informed Neural Networks with Fourier-Enhanced Features

Spectral bias in physics-informed and operator learning: Analysis and mitigation guidelines