1 Introduction↩︎

Partial differential equations (PDEs) govern dynamics in heat transfer, fluid flow, electromagnetics, and other physical systems. Efficient PDE solvers are central to scientific progress. Physics-Informed Neural Networks (PINNs) [@raissi2017physicspart1; @raissi2018numerical] embed PDE constraints into neural network loss functions, enabling both forward and inverse problem solving. PINNs have been applied widely [@cuomo_scientific_2022], but often fail in hard PDE settings, converging to trivial or inaccurate solutions. A key challenge is the imbalance between data-driven and physics-informed loss terms [@wang2021understanding], which causes optimization pathologies. Remedies include curriculum learning [@krishnapriyan_characterizing_2021], adaptive sampling [@daw_mitigating_2023], and spectral/Neural Tangent Kernel (NTK) analyses [@wang2021eigenvector; @wang2022and; @lau2024pinnacle]. Yet these approaches typically succeed only in narrow PDE regimes and lack generality. Meanwhile, the machine learning community has shown that initial weights strongly influence trainability in deep networks [@arpit2019initialize; @schoenholz_deep_2017; @taki2017deep; @yang2017mean]. Standard weight initialization schemes such as Xavier or Kaiming can place networks in chaotic or overly ordered regimes, leading to poor convergence. Despite this, nearly all PINN studies rely on off-the-shelf initializers without task-specific adaptation.

To address this gap, we introduce Learned-Initialization via Gated Layerwise Optimization for PINNs (LIGO-PINN), a framework to systematically initialize PINN network weights that alleviates catastrophic training dynamics. Our contributions are as follows:

• We propose a method to systematically learn the initial PINN weights, to mitigate catastrophic training failures, commonly observed in PINNs, particularly in challenging PDE regimes.

• As part of the learned initialization method, we introduce a novel gated layer-wise optimization (GLO) procedure and show that it produces initial PINN weights that mitigate spectral bias during PINN training.

• Through extensive experiments on three diverse PDE domains, we demonstrate consistent improvements over state-of-the-art baselines, especially in extrapolation settings.

2 Related Work↩︎

We now review major research directions investigating PINN training failures [@cuomo2022scientific].

Temporal Dynamics, Sampling, and Curriculum-Based Methods. A seminal work  [@wang2021understanding] uncovered the imbalance in gradients between data-driven and physics-based loss terms as one of the reasons for catastrophic PINN failures, causing converge to trivial solutions. Early work also highlighted the importance of respecting causal order of PDE dynamics during training, in time-dependent problems [@wang2022respecting]. Building on this, residual-based sampling techniques dynamically reweigh collocation points in regions with high PDE residuals [@daw_mitigating_2023; @toloubidokhti2023dats]. In parallel, curriculum learning approaches schedule training tasks to be ordered by increased difficulty, guiding PINNs from simpler to complex PDE regimes [@krishnapriyan_characterizing_2021]. While these strategies reduce domain errors, they require careful tuning when PDE parameters vary widely and do not generalize well across different PDE systems.

Gradient-Level Interventions Efforts to address the issue of loss imbalance in PINNs have also explored gradient-level solutions. [@kim2021dpm] introduced a Dynamic Pulling Method (DPM) which employs a pseudo-inverse operation to align losses in the same direction, to improve stability, effectiveness of training. [@wang2025gradientalignmentphysicsinformedneural] investigated effect of higher-order optimizers.

Meta-Learning for PINNs. Several recent works have explored meta learning concepts in PDE systems. Psaros et al. [@psaros_meta-learning_2022] uses meta learning to discover optimal loss function for different PINN system. In contrast LIGO-PINN uses meta learning to inform PINN weight initialization. Hyper-LR-PINNs [@cho2023hypernetwork] employ a hyper-network architecture to generate low-rank, task-specific weights for each PDE domain. Although they have trained on similar PDE domains to ours, their experiments have not explored extrapolation regimes as extreme as ours in this paper. [@qin2022metapdelearningsolvepdes] is the closest paper to ours and also employs a model-agnostic meta-learning (MAML) approach for PINN training. However, their evaluation is restricted to simple PDE configurations and unlike our setup, they do not investigate gated optimization strategies or generalization to unseen PDE configurations.

In contrast to other approaches to alleviate catastrophic training failures in PINN, LIGO-PINN is the only approach to systematically explore the effect of learned weight initialization on PINN training dynamics. Further, our gated layer-wise optimization is the first approach of its kind applied in the context of addressing convergence failures in PINNs across diverse PDE settings.

3 Problem Formulation↩︎

We follow the standard PINN setting as in [@raissi2017physicsinformeddeeplearning], also described briefly below; See Appendix 7 for detailed formulation. Let \(\mathbf{z} \in \Omega \subset \mathbb{R}^{d}\) denote the spatial coordinate within domain \(\Omega\) (with boundary \(\partial\Omega\)), \(t \in [0,T]\) time, and \(u(\mathbf{z},t)\) the physical state solving the PDE in Eq. 1 . \[\mathcal{F}\big(u(\mathbf{z},t);\gamma\big) = f(\mathbf{z},t), \quad \mathbf{z} \in \Omega,\;t \in [0,T], \qquad \mathcal{B}\big(u(\mathbf{z},t)\big)=u_b(\mathbf{z},t),\;\mathbf{z} \in \partial\Omega, \label{eq:pde95formulation}\tag{1}\] Here, \(\mathcal{F}\) is the (possibly nonlinear) differential operator representing the underlying physical law, \(\gamma\) are the PDE parameters (e.g., viscosity, diffusivity), and \(f(\mathbf{z},t)\) is a known forcing or source term. The boundary operator \(\mathcal{B}\) enforces boundary or initial conditions, and \(u_b(\mathbf{z},t)\) specifies the prescribed values along \(\partial\Omega\) and the subscript ‘\(b\)’ is added to indicate solution along the boundary \(\partial\Omega\).

A neural network \(\hat{u}_\Theta(\mathbf{z},t)\) parameterized by weights \(\Theta\) is trained to satisfy the PDE at interior (collocation) points and to match boundary/initial data on \(\partial\Omega\). With \(N_{\Omega}\) interior samples \({(\mathbf{z}_i,t_i)}\) and \(N_{\partial\Omega}\) boundary/initial samples \({(\mathbf{z}_j,t_j)}\), the objective is listed in Eq. 2 . \[\label{eq:pinn-total-loss} \mathcal{L}(\Theta) = \lambda_r \frac{1}{N_{\Omega}} \sum_{i=1}^{N_{\Omega}} \Big|\mathcal{F}(\hat{u}_\Theta(\mathbf{z}_i,t_i);\gamma) - f(\mathbf{z}_i,t_i)\Big|^2 + \frac{1}{N_{\partial\Omega}} \sum_{j=1}^{N_{\partial\Omega}} \Big|\hat{u}_\Theta(\mathbf{z}_j,t_j) - u_b(\mathbf{z}_j,t_j)\Big|^2.\tag{2}\] Here, \(\lambda_r\) balances the contributions of the physics residual and boundary losses.

Previous work has tried to alleviate PINN training failures by hyper-parameter tuning of \(\lambda_r\) [@wang2021understanding], learning improved training curricula [@krishnapriyan_characterizing_2021] or by intelligently re-sampling collocation points to focus on hard regions of the domain [@daw_mitigating_2023; @wu2023comprehensive] but has all seen limited success as our investigations of these methods in Sec. 5 reveal.

Complementary to these works, we hypothesize that PINN convergence failures in challenging (a.k.a. hard) PDE domains, can be alleviated, if PINN weights \(\Theta\), instead of being randomly initialized, are systematically initialized prior to training on the hard domain. To this end, we propose a Learned-Initialization mechanism to systematically learn better initial weights, such that PINNs trained with these learned initial weights overcome catastrophic training failures. Our proposed LIGO-PINN method, has two facets (i) Invariance Encoding (ii) Gated Layer-wise Optimization.

Invariance Encoding (IE). If the goal of the PINN is to learn a challenging PDE domain governed by PDE parameter \(\gamma_{\mathrm{hard}}\), IE seeks to first train a randomly initialized PINN model on a set of relatively easier configurations of the PDE with parameters \(\Gamma_{\mathrm{easy}} = \{\gamma_{1},\dots,\gamma_{K}\}\). This enables a distillation of the invariant physics across \(\Gamma_{\mathrm{easy}}\) into the PINN weights \(\Theta\). Our approach for IE is grounded in foundational works of meta‐learning [@finn_model-agnostic_2017]. IE begins with randomly initialized global PINN weights \(\Theta_0\) and a pool of \(K\) tasks in \(\Gamma_{\mathrm{easy}}\). The IE learning process is carried out for \(\mathcal{R}\) iterations \(\{1,2,\dots,\mathcal{R}\}\) and at each training iteration \(r\), tasks \(\{\mathcal{T}_1,\dots,\mathcal{T}_k\}\) are randomly sampled from the task distribution \(p(\Gamma_{\mathrm{easy}})\) of easy tasks and \(k\) identical copies \(\{\Theta'_{r,1},\dots,\Theta'_{r,k}\}\) of the current global PINN weights \(\Theta_r\) are created. Each \(\Theta'_{r,i}\) is then optimized via. gradient descent, employing data from task \(\mathcal{T}_i\) at iteration \(r\) for ‘J’ inner iterations (for simplicity we assume ‘J’ = 1). Let \(\Theta''_{r,i}\) represent the updated version of \(\Theta'_{r,i}\) after ‘J’ inner gradient-descent update iterations. The common representations across all \(k\) task models \(\{\Theta''_{r,1},\dots,\Theta''_{r,k}\}\) is distilled into \(\Theta_r\) via. a ‘meta-update’ as indicated in Eq. [eq:meta95update] where \(\mathcal{L}_{\mathcal{T}_i}(\Theta''_{r,i})\) represents the loss in Eq. 2 calculated w.r.t parameters \(\Theta''_{r,i}\) with data from task \(\mathcal{T}_i\) and \(\eta_{\mathrm{outer}}\) represents the learning rate.

\[\begin{align} \Theta_{r+1} = \Theta_{r} - \eta_{\mathrm{outer}}\nabla_{\Theta_r}\mathcal{L}_{\mathrm{meta}};\qquad\qquad \nabla_{\Theta_r}\mathcal{L}_{\mathrm{meta}} =\nabla_{\Theta_r}\sum_{i=1}^k\mathcal{L}_{\mathcal{T}_i}(\Theta''_{r,i})\label{irqyahus} \end{align}\tag{3}\]

Gated Layer-wise Optimization (GLO). Traditional meta-learning approaches [@finn_model-agnostic_2017; @rajeswaran2019meta] are not physics-informed and treat all tasks in \(\Gamma_{\mathrm{easy}}\) as resulting in similar training dynamics across all layers. Thus, the meta-update in Eq. [eq:meta95update] is a simple summation of all task gradients. However, based on previous work [@wang2021understanding], we know that non-trivial dynamics between the physics-based loss and data-driven loss affects different layers in the neural network differently. To address these gaps that manifest uniquely in the physics-informed PDE modeling context, we augment the IE procedure with the novel GLO mechanism by introducing a gated update of layers during each meta-update, such that a gating mechanism \(g(r)\), unlocks deeper layers, exposing their parameters to gradient-based optimization, gradually, over the training iterations from 1 to \(\mathcal{R}\).

Let us consider that the global PINN is an ‘L’ layer neural network whose parameters at iteration \(r\) are represented as \(\Theta_r = \{\theta_r^0,\dots,\theta_r^{L-1}\}\) with \(\theta^i_r\) representing the weights corresponding to layer \(i\) at iteration \(r\). GLO allows us to prioritize learning shallow layers prior to deeper layers providing fine-grained (layer level) control over training dynamics in contrast to traditional IE approaches. Specifically, at each iteration \(r \in \{1,2,\dots,\mathcal{R}\}\), the (binary) gating value for layer \(l\) is calculated with a gating mechanism detailed in Eq. 4 . \[g^l(r) = \begin{cases} 1, & l \;\le\; \bigl\lfloor\tfrac{r\,L}{\mathcal{R}}\bigr\rfloor,\\ 0, & \text{otherwise}, \end{cases} \quad l=0,\dots,L-1. \label{eq:linear95gating}\tag{4}\] The gating schedule begins with the first layer unmasked, followed by linear unmasking of deeper layers as training progresses, ensuring localized (and more stable) effects of non-trivial loss dynamics.

Although the proposed \(g(r) \in \{0,1\}^L\), results in a linear gating schedule, GLO can easily admit any general gating mechanism.

\[\theta^{l}_{r+1} =\theta^{l}_r -\eta_{\mathrm{outer}}\,g^l(r)\,\nabla_{\theta_r^l}\mathcal{L}_{\mathrm{meta}}, \quad l=0,\dots,L-1, \label{eq:glo95layer95wise95meta95update}\tag{5}\]

If \(g(r) = [g^0(r),\,\dots,\,g^{L-1}(r)]\) and each meta-update gradient vector is represented as \(\nabla_{\Theta_r}\mathcal{L}_{\mathrm{meta}} =[\nabla_{\theta^0_r}\mathcal{L}_{\mathrm{meta}},\,\dots,\,\nabla_{\theta_r^{L-1}}\mathcal{L}_{\mathrm{meta}}]\), then the meta-update at iteration \(r+1\) for layer \(l\) is detailed in Eq. 5 . Here, each \(g^l(r)\) controls the meta-update strength of layer \(l\).

Alg. 1 details the entire procedure of our proposed LIGO-PINN framework. The learned weights \(\Theta_\mathcal{R}\) resulting from Alg. 1, can be employed as initial weights for PINN training on \(\gamma_{\mathrm{hard}}\) (i.e., the hard PDE context).

4 Dataset Description and Experimental Setup↩︎

We evaluate our proposed LIGO-PINN model on three diverse PDE domains, described below.

1D Convection. The 1D convection PDE (Eq. 6 ) models transport of a scalar field \(u(z,t)\) with constant velocity. In Eq. 6 , \(z\in[0,2\pi]\) denotes the spatial coordinate and \(t \in[0,1]\) denotes the temporal coordinate. \(\beta\) is the advection coefficient, \(u_t\) and \(u_z\) denote partial derivatives in time and space. \(u(z,0)=sin(z)\) defines the initial condition and periodic boundaries impose \(u(0,t)=u(2\pi,t)\). \[\mathcal{N}[u](z,t) := u_t(z,t) + \beta\,u_z(z,t) = 0, \label{eq:conv}\tag{6}\]

This problem represents convection-dominated motion, where high-speed phase propagation and increasingly oscillatory space-time structure make PINN training challenging. In line with prior work [@krishnapriyan_characterizing_2021], we consider \(\Gamma_{\mathrm{easy}}\) as \(\beta\in\{5,10,15,20,25\}\) and hard regime as \(\beta\in\{40,50,60,70,80\}\). Selection methodology for easy, hard tasks is detailed in Appendix 8.2.

2D Helmholtz. The steady-state Helmholtz PDE (Eq. 7 ) describes oscillatory wave fields on a 2D domain. Let \(z=(x_1,x_2)\) and let \(\Delta=\partial_{x_1x_1}+\partial_{x_2x_2}\) denote the Laplace operator. Eq. 7 is defined on the domain \([-1,1]^2\) with Dirichlet boundary condition \(u=0\) on \(\partial\Omega\). Using the analytic solution \(u(z)=\sin(a_1\pi x_1)\sin(a_2\pi x_2)\), the source term becomes \(q(z)=\big[k^2-(a_1\pi)^2-(a_2\pi)^2\big]\,u(z)\). \[\mathcal{N}[u](z) := \Delta u(z) + k^2 u(z) - q(z) = 0, \label{eq:helm}\tag{7}\] The coefficients \((a_1,a_2)\) control the oscillation frequency along x and y direction respectively, and larger values induce more frequent oscillations posing greater challenge for PINNs. \(\Gamma_{\mathrm{easy}}\) is considered \(\{(1,1),(1,2),(1,3),(1,4),(1,5)\}\) based on the case considered in [@wang2021understanding] and hard regime is selected to have significantly more challenging propagation dynamics and both isotropic and anisotropic propagation dynamics are tested in the hard regimes \(\{(4,4),(4,5),(5,5),(4,6),(5,6),(6,6)\}\). Our work has significantly expanded the investigation regime, compared to previous PINN-based investigations of 2D Helmholtz. Selection methodology for easy, hard tasks is detailed in Appendix 8.2..

2D Incompressible Navier–Stokes. We finally consider modeling the challenging 2D viscous flow over a cylinder, using the non-dimensional 2D incompressible Navier–Stokes equations, represented in Eq. 8 in compact momentum-continuity form. Here \(\mathbf{u}(x,y,t)\) is the velocity field, \(p(x,y,t)\) is the scalar pressure field, \(\nu\) is the kinematic viscosity (parameterized through the Reynolds number). Domain \((x,y,t)\in\Omega\times[0,T]\), where \(\Omega\) contains a unit-diameter cylinder. \[\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \nabla p - \nu \nabla^{2}\mathbf{u} = 0, \qquad \nabla\cdot\mathbf{u}=0.\label{eq:ns}\tag{8}\] This system exhibits vortex shedding and transitional flow behavior, and the coupling between momentum and incompressibility creates multiple interacting residual terms, making it one of the most challenging PDE families for PINNs. The complexity of the dynamics increases with increase in Reynolds number (\(Re\)). We consider \(\Gamma_{\mathrm{easy}}\) to be \(\mathrm{Re}\in\{100,200,300,400,500\}\) where relatively more steady, laminar dynamics are exhibited and in line with [@lee2025physics], we consider \(\mathrm{Re}\in\{600,800,1000\}\) to be hard regime where less steady dynamics are exhibited.

To ensure a fair evaluation, all PINN baseline methods are trained for 56K epochs on the target domain \(\Gamma_{\mathrm{hard}}\), while LIGO-PINN is trained only for 50K epochs on the target domain. This difference accounts for the 6K iterations employed for the learned initialization procedure which is absent in baseline methods which are randomly initialized. All other training parameters (detailed in Appendix 8 for full reproducibility) are identical including model architecture, training data, optimization settings.

5 Results & Discussion↩︎

In this section, we demonstrate with rigorous quantitative and qualitative experiments, the effects of catastrophic failures on PINN performance and how these failures can be alleviated via Learned Initialization. To validate our approach, we compare LIGO-PINN against a suite of state-of-the-art models, that have been designed specifically to overcome catastrophic failures of PINNs in challenging PDE domains; the full descriptions of these baselines can be found in Appendix 9.2. Specifically, we undertake this investigation by answering three research questions (RQ):

RQ1. (Qualitative) How does learned initialization alleviate catastrophic PINN training dynamics?

RQ2. (Quantitative) How does LIGO-PINN compare with other state-of-the-art (SoTA) approaches proposed to alleviate PINN catastrophic failures?

RQ3. (Ablation) What is the effect of invariance encoding and gated layer-wise optimization in LIGO-PINN?

5.1 (RQ1) How does learned initialization alleviate catastrophic PINN training dynamics?↩︎

We answer this question in three stages. First, we demonstrate qualitatively the catastrophic failures PINNs undergo in challenging PDE contexts. Second, we investigate how the training losses of the PINNs evolve in challenging domains. Third, we characterize the effect of the catastrophic training dynamics by analyzing the loss landscape neighborhood of the converged models. In each case, we contrast the analysis of a randomly initialized PINN, with LIGO-PINN behavior in the same context.

Qualitative Failure Modes. In Fig. 2a, we demonstrate the solution predicted by a PINN (with random weight initialization) when trained on a challenging 1D convection domain with \(\beta = 70\). It can be clearly seen that the PINN model (center plot in Fig. 2a) has failed to estimate a good solution of the domain. The perfect solution of the domain can be seen in the left plot (labeled ‘Analytical Solution’) of Fig. 2a. We conduct a similar qualitative investigation of a separate PDE domain 2D-Helmholtz comprising two PDE parameters \(a_1, a_2\), with \((a_1 = 6, a_2 = 6)\) in Fig. 2b. Once again, we notice that the PINN model (center plot of Fig. 2b) is unable to provide a good estimate of the challenging PDE domain solution (i.e., ‘Analytical Solution’ represented on the left plot in Fig. 2b). Turning our attention to the performance of LIGO-PINN in each of these two contexts (i.e., right plot in each of Fig. 2a and Fig. 2b), we notice that the estimation of LIGO-PINN closely aligns with the analytical solution of the corresponding PDE domain. We now highlight that the only difference between the PINN and LIGO-PINN procedures is the weight initialization. Specifically, PINN weights have been initialized randomly (i.e., Xavier initialization) while the weights of LIGO-PINN have been initialized by our proposed learned initialization method. Thus, we infer that the learned initialization strategy in LIGO-PINN successfully alleviates catastrophic failures to which traditional PINNs are susceptible, when modeling complex PDE dynamics.

Extension to Multi-Physics. We also investigate LIGO-PINN performance on 2D incompressible, Navier-Stokes cylinder wake at high Reynolds numbers. A challenging multi-physics setting. We investigate both forward (pressure prediction) and inverse settings (recovering viscous and advection coefficients of Navier-Stokes equations).

Consistent with our findings on 1D-Convection and 2D-Helmholtz, we see from Fig. 3, that LIGO-PINN yields substantially lower errors and qualitatively faithful pressure fields (Fig. 3 (f)) compared to randomly initialized PINN (Fig. 3 (e)). Ground truth pressure field is depicted in Fig. 3 (d). This highlights the advantage of learned initialization to model complex fluid dynamics.

Training Dynamics Analysis. We complement the qualitative analysis in Fig. 2 with an investigation (Fig.4) of the training dynamics in the same PDE contexts. Specifically, we highlight in Fig.4, the evolution of the total training loss (Eq. 2 ) of PINNs and LIGO-PINN on the same 1D convection, 2D Helmholtz and 2D Navier-Stokes domains. We note that the training loss plots in Fig.4 has been averaged over five independent training runs to capture the expected training dynamics.

For the 1D convection problem (\(\beta = 70\)) the result depicted in Fig. 4a, shows that random initialization of the PINN leads to an early plateau in the loss curve, whereas our learned initialization method (LIGO-PINN) achieves a steady, monotonic decrease, signaling more effective exploration of the solution manifold and ultimately, convergence to a loss value an order of magnitude lower than the corresponding randomly initialized PINN solution. A similar pattern emerges on the 2D Helmholtz and 2D Navier-Stokes problem \((a_1=6,\,a_2=6)\) as shown in Fig. 4b,c. We see that the randomly initialized PINN exhibits premature convergence to a sub-optimal solution while, once again LIGO-PINN achieves a loss several orders of magnitude lower. This demonstrates the ability of our learned initialization based solution (LIGO-PINN) to achieve superior convergence dynamics by overcoming loss surface barriers [@fort2020deep] and hence avoiding catastrophic failures.

Loss Landscape Analysis. Although the loss curves in Fig. 4 show a clear contrast between random and learned initialization, they do not fully reveal how training dynamics affect parameter optimization. To clarify this, we examine the loss landscape neighborhood of each converged model using a Hessian‐informed perturbation analysis [@krishnapriyan_characterizing_2021; @B_ttcher_2024]. Fig. 5 presents results for the 1D convection problem (2D Helmholtz in Appendix 9.2). Under random initialization (Fig. 5a), the converged PINN exhibits a highly rugged surface with deep valleys and sharp ridges, driven mainly by the residual term (Fig. 5a.iii). The presence of such ravines is known to hinder convergence to strong local minima [@wu2017towards; @hochreiter1997flat]. In contrast, the data loss landscape is smoother but still challenging to optimize. With LIGO-PINN, the converged landscape is markedly smoother without the sharp peaks of Fig. 5a. Fig. 5b.i shows the total loss surface, while Figs. 5b.ii–b.iii show data and residual components, all with smoother geometry compared to their counterparts resulting from a randomly initialized PINN. Viewed together with Fig. 4a, these results illustrate that random initialization produces residual-dominated peaks that trap the optimizer in narrow basins, whereas LIGO-PINN yields a smoother landscape that supports improved convergence to better local minima.

5.2 (RQ2) How do LIGO-PINNs compare with other state-of-the-art (SoTA) approaches proposed to alleviate PINN catastrophic failures?↩︎

Recent efforts have begun to investigate catastrophic failures in PINNs, and in this section we compare our learned-initialization-based LIGO-PINN model against six comparative PINN baselines, including vanilla PINNs and recent methods proposed to improve PINN training stability [@krishnapriyan_characterizing_2021; @wang2021understanding; @daw_mitigating_2023; @qin2022metapdelearningsolvepdes; @cho2023hypernetwork]. We provide quantitative comparisons using the mean-absolute error (MAE) metric, across three benchmark PDE systems: the 1D Convection equation, the 2D Helmholtz equation, and the 2D incompressible time-dependent Navier-Stokes system [@raissi2017physicspart1]. In all cases, we evaluate models under challenging PDE settings, and extend our analysis regime far beyond what existing papers have investigated. For Navier-Stokes, we report both forward pressure reconstruction and inverse coefficient estimation in the main text.

1D Convection System Table 1 shows that LIGO-PINN achieves the lowest average MAE across the extrapolation settings and obtains the best result in four out of five tested regimes.  LIGO-PINN achieves the strongest average performance, reducing average MAE by 96.5% relative to vanilla PINNs and by 68.9% relative to the strongest non-ours baseline, PINN-Dynamic. Across all six comparative baselines, the MAE reduction ranges from 68.9% to 96.5%.

2D Helmholtz System. The 2D Helmholtz system presents a more challenging oscillatory and multi-parameter setting. As shown in Table 2, the baseline methods remain substantially less accurate than LIGO-PINN across all tested domain parameters. On average, LIGO-PINN reduces MAE by 94.3% relative to vanilla PINNs and by 92.6% relative to the strongest non-ours baseline, PINN-Dynamic. Across all six comparative baselines, the MAE reduction ranges from 92.6% to 95.5%. These results indicate that learned initialization is particularly effective for high-frequency PDE regimes where standard PINN training and existing mitigation strategies remain unstable.

5pt

2D Incompressible Navier–Stokes. Table 3 reports the pressure prediction results across Reynolds numbers. In this fluid dynamics setting, LIGO-PINN achieves an average MAE of \(0.0096\), compared with \(0.0517\) of PINN (dynamic), (the strongest baseline), and \(0.0791\) for vanilla PINNs (i.e., PINN (fixed)). This corresponds to an error reduction of 81.4% relative to the strongest baseline and 87.9% relative to vanilla PINNs. Across all six comparative baselines, the MAE reduction ranges from 81.4% to 97.4%. These results show that learned initialization remains effective even in complex time-dependent fluid dynamics problems.

In addition to the forward pressure reconstruction setting, we further evaluate all methods on the inverse formulation of the 2D incompressible Navier–Stokes system. In this setting, the advection and viscosity coefficients \((\lambda_{\mathrm{adv}}, \lambda_{\mathrm{visc}})\) are treated as unknown parameters and jointly inferred during training from sparse observations. This inverse problem is particularly challenging due to the strong coupling between nonlinear momentum dynamics and incompressibility constraints, especially at higher Reynolds numbers. Table 4 reports the coefficient estimation errors across Reynolds numbers \(\mathrm{Re}\in\{600,800,1000\}\). We observe that standard PINN baselines exhibit unstable parameter recovery and large estimation errors as the flow dynamics become increasingly complex. In contrast, LIGO-PINN consistently achieves lower coefficient estimation error across all evaluated regimes, demonstrating that learned initialization improves not only forward field reconstruction but also inverse parameter identification in challenging multi-physics PDE systems. For the inverse problem involving estimation of the advection (\(\lambda_{\mathrm{adv}}\)) and viscosity (\(\lambda_{\mathrm{visc}}\)) coefficients, LIGO-PINN yields an 80.2%, and 45.9% performance improvement over the strongest baseline (HyperLR in the inverse modeling case).

3pt

5.3 (RQ3) What is the effect of invariance encoding and gated layer-wise optimization in LIGO-PINN? (Ablation Analysis)↩︎

We analyze the effect of Invariance Encoding (IE) and Gated Layer-wise Optimization (GLO) through ablation variants: ‘LIGO-PINN w/o (IE,GLO)’ (random initialization) and ‘LIGO-PINN w/o GLO’ (MAML-style initialization). We focus on the challenging 2D Navier-Stokes setting, a multi-physics, highly coupled domain where PINN failures are most pronounced. The ablation analysis is conducted in two stages, (i) Empirical Evaluation (MAE) (ii) Spectral Evaluation of Absolute Error Fields.

Empirical Evaluation. Fig. 6 (a) shows that neither simple random initialization (‘w/o (IE,GLO)’) nor just incorporating IE (‘w/o GLO’) improves performance on its own. Effective performance is only demonstrated by LIGO-PINN where both IE and GLO simultaneously influence learning initialization. This highlights the critical role of gated, layer-wise parameter optimization and invariance encoding in learning initial weights that enable stable and accurate convergence in challenging PDE domains.

Spectral Evaluation. Fig. 6 (b) shows spectral density of the absolute error fields of each ablation variant (lower is better). LIGO-PINN consistently achieves lower error across both low and high frequencies, while ablation variants exhibit pronounced high-frequency error (highlighting spectral bias). Specifically, this indicates improved capture of fine-scale structures in fluid dynamics by LIGO-PINN compared to variants without IE or GLO. Results on 1D Convection and 2D Helmholtz exhibit similar trends and detailed in Appendix 9.6.

Additional Results. We also successfully verify that LIGO-PINN performance translates to complex unstructured 3D domains (results in Appendix 9.1). For completeness, we also report computational time and memory footprint details in Appendix 8.5. Our results depict that LIGO-PINN scales comparably to standard PINNs aside from a one-time meta-training overhead. Appendix 9.7 presents sensitivity analyses, to the number and range of sampled tasks during learned initialization phase.

6 Conclusion↩︎

This work shows that learned initialization can serve as an effective mechanism for reducing catastrophic failure in PINNs. Rather than relying only on task-specific training heuristics such as curriculum design, adaptive sampling, or loss reweighting, LIGO-PINN learns a weight initialization that places the PINN in a more favorable optimization regime before training begins, leading to more stable training and substantially lower errors on challenging extrapolation problems. Across the evaluated PDE systems, LIGO-PINN reduces average error by 91.5% over six comparative PINN baselines and improves over the strongest baseline in each domain by an average of 81%. The ablation results further indicate that both invariance encoding and gated layer-wise optimization are important contributors to this behavior, reducing volatility and improving spectral stability. These findings suggest that learned initialization is a promising direction for more robust PINN training, especially in regimes where standard training often breaks down. Finally, the results on 3D unstructured domains indicate that the benefits of LIGO-PINN are not ONLY limited to structured 2D settings and generalize to complex domains.

Limitations and Future Work. While LIGO-PINN demonstrates strong performance across several challenging PDE systems, the proposed invariance encoding stage requires access to related PDE configurations during initialization. In addition, although the empirical analyses strongly suggest improved optimization behavior and reduced spectral bias, a deeper theoretical understanding of the underlying mechanisms will strengthen the framework and will be explored in the future. Finally, extending LIGO-PINN to real-world large-scale multi-physics and industrial simulation settings is a promising avenue for future investigation.

7 Detailed formulation of PINN↩︎

In this section, we first present the general formulation of Physics-Informed Neural Networks (PINNs), and then demonstrate its instantiation on two canonical PDEs: the 1D Convection Equation and the 2D Helmholtz Equation.

7.1 Physics Informed Neural Network↩︎

A general form of a PDE can be expressed as: \[\mathcal{F}(u(\mathbf{z}); \mathbb{\gamma}) = f(\mathbf{z}), \quad \mathbf{z} \in \Omega, \label{eq:physicspde}\tag{9}\] \[\mathcal{B}(u(\mathbf{z})) = g(\mathbf{z}), \quad \mathbf{z} \in \partial \Omega.\]

Here:

• \(\Omega \subset \mathbb{R}^d\) represents the spatial domain, while \(\partial \Omega\) denotes its boundary.

• \(\mathbf{z} = [x_1, x_2, \dots, x_d, t]^\top \in \mathbb{R}^{d+1}\) encapsulates space-time coordinates.

• \(u\) is the unknown solution we aim to approximate.

• \(\mathbb{\gamma}\) denotes physical parameters of the system.

• \(f(\mathbf{z})\) encodes system-specific data or forcing terms.

• \(\mathcal{F}\) is a (potentially nonlinear) differential operator characterizing the physical laws of the system.

• \(\mathcal{B}\) is a boundary operator that enforces initial and boundary conditions, which can be of Dirichlet, Neumann, or periodic types.

The PINN framework leverages automatic differentiation to embed the residuals of \(\mathcal{F}\) and \(\mathcal{B}\) into a composite loss function. This loss function penalizes deviations from the governing equations and boundary conditions, ensuring that the solution adheres to the physical laws while simultaneously fitting observed data.

Neural networks are powerful tools due to their universal approximation property, which allows them to approximate any continuous function given sufficient capacity. In the context of PINNs, we aim to approximate the solution \(u(\mathbf{z})\) of the governing equations using a neural network \(\hat{u}_\theta(\mathbf{z})\), where \(\theta\) denotes the learnable parameters of the neural network. The approximation can be written as: \[\hat{u}_\theta(\mathbf{z}) \approx u(\mathbf{z}),\] where \(\mathbf{z} \in \mathbb{R}^{d+1}\) represents the space-time coordinates.

To achieve this, the PINN framework constructs a composite loss function consisting of:

1. PDE Residual Loss: This enforces the governing physical laws as specified by the PDE.

2. Boundary/Data Loss: This ensures compliance with boundary or initial conditions, as well as data consistency when available.

7.1.0.1 PDE Residual Loss

Using automatic differentiation tools, the neural network computes derivatives efficiently with respect to its inputs, enabling evaluation of the PDE residual. The residual is given by: \[\mathcal{R}(\mathbf{z}) = \mathcal{F}(\hat{u}_\theta(\mathbf{z}); \mathbb{\gamma}) - f(\mathbf{z}),\] where \(\mathcal{F}\) represents the differential operator, \(f(\mathbf{z})\) is the system-specific forcing term, and \(\mathbb{\gamma}\) are the physical parameters. The residual loss over the domain \(\Omega\) is formulated as: \[\mathcal{L}_{{r}} = \frac{1}{N_\Omega} \sum_{i=1}^{N_\Omega} \left| \mathcal{R}(\mathbf{z}_i) \right|^2,\] where \(N_\Omega\) represents the number of collocation points in the domain, and \(\mathbf{z}_i \in \Omega\) are the sampled points.

7.1.0.2 Boundary/Data Loss

For enforcing boundary or initial conditions, we define a loss term that penalizes deviations from the specified boundary values \(g(\mathbf{z})\). For Dirichlet conditions, this can be written as: \[\mathcal{L}_{{b}} = \frac{1}{N_{\partial \Omega}} \sum_{j=1}^{N_{\partial \Omega}} \left| \mathcal{B}(\hat{u}_\theta(\mathbf{z}_j)) - g(\mathbf{z}_j) \right|^2,\] where \(N_{\partial \Omega}\) represents the number of points on the boundary \(\partial \Omega\), and \(\mathbf{z}_j \in \partial \Omega\) are the boundary points.

7.1.0.3 Total Loss

The total loss for training the PINN combines the residual loss and the boundary/data loss, weighted by respective coefficients \(\lambda_{{r}}\) and \(\lambda_{{d}}\): \[\mathcal{L}_{{total}} = \lambda_{{r}} \mathcal{L}_{{r}} + \lambda_{{d}} \mathcal{L}_{{b}}. \label{eq:pinn95total95loss}\tag{10}\]

We calculate the total loss by adding the residual and data losses. Here, \(\lambda_{{r}}\) and \(\lambda_{{d}}\) are hyperparameters that control the relative contribution of the residual and data losses during training. The residual loss acts as a regularization term, ensuring that the learned solution adheres to physically consistent laws.

To train the neural network, we solve the optimization problem: \[\theta^* = \arg\min_{\theta} \mathcal{L}_{{total}}. \label{eq:optim95loss}\tag{11}\]

We aim to formulate a Physics-Informed Neural Network (PINN) by studying two fundamental classification of systems in physics that are crucial for understanding and modeling the physical world. These systems are categorized based on whether their properties depend on time:

Time-Independent Systems These systems are used to study the steady-state behavior of a system, where properties such as energy, momentum, or other conserved quantities remain constant over time. Time-independent systems help us analyze static or equilibrium conditions, providing insights into the inherent properties of the system without considering temporal evolution.

Time-Dependent Systems These systems are used to study how the state of a system evolves over time. By understanding the temporal dynamics, we can gain insights into transient phenomena, system responses, and changes in state variables due to external forces or intrinsic behaviors.

To explore these concepts, we explore the following two systems.

7.2 1D Convection Equation↩︎

This system represents a classic example of a time-dependent system. We consider a one-dimensional convection system to demonstrate the application of Physics-Informed Neural Networks (PINNs). The system models the transport of heat or mass in a medium, where the primary variable \(u(x, t)\) evolves over space and time. The spatial domain is \(x \in [0, 2\pi]\), and the time domain is \(t \in [0, 1]\). This setup allows for studying convection behavior under periodic boundary conditions, often encountered in cyclic systems. This formulation is adapted from  [@krishnapriyan_characterizing_2021]

The governing partial differential equation (PDE) is given by: \[\frac{\partial u}{\partial t} + \beta \frac{\partial u}{\partial x} = 0, \quad x \in \Omega, \, t \in [0, T],\] with the initial condition: \[u(x, 0) = h(x), \quad x \in \Omega.\]

Here:

• \(\beta\) is the convection coefficient (physical parameter of system \(\gamma\) in 9 ),

• \(h(x)\) specifies the initial condition for \(u(x, t)\),

• \(\Omega\) denotes the spatial domain.

To evaluate the performance of PINNs, we generate ground truth data using an analytical solution. For constant \(\beta\) and periodic boundary conditions, the analytical solution is: \[u_{{analytical}}(x, t) = \mathcal{F}^{-1} \left[ \mathcal{F}(h(x)) e^{-i\beta k t} \right],\] where:

• \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) are the Fourier transform and its inverse, respectively,

• \(i = \sqrt{-1}\) is the imaginary unit,

• \(k\) represents the frequency in the Fourier domain.

For this problem, we assume: \[h(x) = \sin(x), \quad u(0, t) = u(2\pi, t),\] implying periodic boundary conditions.

7.2.0.1 Residual Loss

The residual loss for this PDE is given by: \[\mathcal{L}_{{r}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left( \frac{\partial \hat{u}}{\partial t} + \beta \frac{\partial \hat{u}}{\partial x} \right)^2,\] where:

• \(N_f\): Number of collocation points sampled from the interior of the domain \(\Omega \times [0, T]\)

• \(\hat{u}\) is the predicted solution,

• \(\frac{\partial \hat{u}}{\partial t}\): Temporal derivative of the predicted solution,

• \(\frac{\partial \hat{u}}{\partial x}\): Spatial derivative of the predicted solution,

• \(\beta\):Convection coefficient which is the parameter describing the physical system

7.2.0.2 Data Loss

The data loss consists of two components: the loss enforcing the initial condition and the loss enforcing the boundary conditions. It is defined as: \[\mathcal{L}_{{d}} = \mathcal{L}_{{ic}} + \mathcal{L}_{{bc}},\] where: \[\mathcal{L}_{{ic}} = \frac{1}{N_u} \sum_{i=1}^{N_u} \left( \hat{u}(x_i, 0) - g(x_i, 0) \right)^2,\] \[\mathcal{L}_{{bc}} = \frac{1}{N_b} \sum_{j=1}^{N_b} \left( \hat{u}(0, t_j) - \hat{u}(2\pi, t_j) \right)^2.\]

Here:

• \(\mathcal{L}_{{ic}}\): Initial condition loss, ensuring the predicted solution \(\hat{u}(x, 0)\) matches the data \(g(x, 0) = h(x)\),

• \(\mathcal{L}_{{bc}}\): Boundary condition loss, ensuring periodic boundary conditions \(\hat{u}(0, t) = \hat{u}(2\pi, t)\),

• \(N_u\): Number of points sampled from the initial condition,

• \(N_b\): Number of points sampled from the boundary condition,

• \(g(x_i, 0) = h(x_i)\): Known solution at initial condition points,

• \(\hat{u}(0, t_j)\) and \(\hat{u}(2\pi, t_j)\): Predicted solutions at boundary points.

7.2.0.3 Total Loss

The total loss combines the above components: \[\mathcal{L}(\theta) = \lambda_{{r}} \mathcal{L}_{{r}} + \lambda_{{d}} \mathcal{L}_{{d}},\] where:

• \(\lambda_{{r}}, \lambda_{{d}}\) are weighting coefficients balancing the contributions of the respective terms.

7.3 2D Helmholtz Equation↩︎

We study the vibration states of a two-dimensional membrane using the 2D Helmholtz equation. This equation is widely used to describe stationary wave fields, such as vibrations in membranes, sound waves, and electromagnetic waves. Specifically, we consider a square membrane with spatial coordinates \(x_1, x_2 \in [-1, 1]\) and analyze its steady-state behavior under fixed boundary conditions. The Helmholtz equation helps identify resonance patterns and steady-state vibration modes in such systems. The formulation is adapted from [@wang2021understanding].

The governing PDE for the 2D Helmholtz equation is: \[\frac{\partial^2 u}{\partial x_1^2} + \frac{\partial^2 u}{\partial x_2^2} + k^2 u = q(x_1, x_2), \quad (x_1, x_2) \in \Omega,\] where \(\Omega = [-1, 1] \times [-1, 1]\) is the spatial domain.

The boundary of the domain, \(\partial \Omega\), represents the fixed edges of the membrane. Consequently, the displacement \(u(x_1, x_2)\) at the boundaries is zero, giving rise to the Dirichlet boundary condition: \[u(x_1, x_2) = 0, \quad (x_1, x_2) \in \partial \Omega.\]

We specifically choose \(k^2 \neq a_1^2 + a_2^2\) to ensure that the selected \(k\) does not correspond to an eigenmode of the domain, resulting in a non-zero source term \(q(x_1,x_2)\). Thus, the PDE system is inherently inhomogeneous.

Here:

• \(k\) is the wavenumber and is kept fixed at \(k=1\).

• \(a_1, a_2 \in \mathbb{Z}^+\) are integer mode numbers representing the spatial oscillation patterns. These serve as the primary PDE parameters affecting solution complexity and difficulty.

• \(q(x_1, x_2)\) is the resulting source term that characterizes external excitation or forcing within the membrane, driving the system response.

For this setup, we use an analytical solution of the form: \[u(x_1, x_2) = \sin(a_1 \pi x_1) \sin(a_2 \pi x_2),\] which satisfies the zero-displacement boundary condition. The corresponding source term is derived as: \[q(x_1, x_2) = \left[ k^2 - (a_1 \pi)^2 - (a_2 \pi)^2 \right] \sin(a_1 \pi x_1) \sin(a_2 \pi x_2).\]

7.3.0.1 Residual Loss

The residual loss quantifies deviations from the Helmholtz equation at collocation points in \(\Omega\): \[\mathcal{L}_{\text{residual}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left( \frac{\partial^2 \hat{u}}{\partial x_1^2} + \frac{\partial^2 \hat{u}}{\partial x_2^2} + k^2 \hat{u} - q(x_{1,i}, x_{2,i}) \right)^2,\] where:

• \(N_f\): Number of collocation points sampled from the interior of the domain \(\Omega\),

• \(\hat{u}\): Predicted solution obtained from the neural network,

• \(q(x_{1,i}, x_{2,i})\): Source term evaluated at the collocation point \((x_{1,i}, x_{2,i})\).

7.3.0.2 Boundary Condition Loss

The boundary loss enforces zero displacement at the edges of the membrane, consistent with the fixed boundary condition: \[\mathcal{L}_{\text{boundary}} = \frac{1}{N_b} \sum_{i=1}^{N_b} \hat{u}(x_{1,i}, x_{2,i})^2,\] where \(N_b\) is the number of points sampled from the boundary \(\partial \Omega\).

7.3.0.3 Total Loss

The total loss combines the residual and boundary condition losses: \[\mathcal{L}(\theta) = \lambda_{r} \mathcal{L}_{\text{residual}} + \lambda_{b} \mathcal{L}_{\text{boundary}},\] where \(\lambda_{r}\) and \(\lambda_{b}\) are weighting coefficients balancing the contributions of the residual and boundary losses.

7.4 2D Incompressible Navier-Stokes (with learnable coefficients)↩︎

We consider the two-dimensional incompressible Navier-Stokes equations in non-dimensional form, augmented with learnable coefficients for the advection and viscosity terms: \[\begin{align} u_t + \lambda_{\text{adv}} (u u_x + v u_y) &= -\,p_x + \lambda_{\text{vis}} (u_{xx} + u_{yy}), \tag{12}\\ v_t + \lambda_{\text{adv}} (u v_x + v v_y) &= -\,p_y + \lambda_{\text{vis}} (v_{xx} + v_{yy}), \tag{13}\\ u_x + v_y &= 0. \tag{14} \end{align}\] Here, \(u(t,x,y)\) and \(v(t,x,y)\) are the velocity components, \(p(t,x,y)\) is the pressure, and \(\lambda_{\text{adv}}, \lambda_{\text{vis}}\) are trainable scalars. In a forward setting, \(p\) is predicted by the network but not supervised; in an inverse setting, \(\lambda_{\text{adv}}\) and \(\lambda_{\text{vis}}\) are learned from data.

7.4.0.1 Network Outputs

Let a neural network \(\mathcal{N}_\theta(t,x,y)\) predict \[(\hat{u}, \hat{v}, \hat{p}) = \mathcal{N}_\theta(t,x,y),\]

7.4.0.2 Residual Loss

Define the PDE residuals at collocation points \(\{(t_i,x_i,y_i)\}_{i=1}^{N_f}\) by \[\begin{align} f_u &= \hat{u}_t + \lambda_{\text{adv}} (\hat{u}\,\hat{u}_x + \hat{v}\,\hat{u}_y) + \hat{p}_x - \lambda_{\text{vis}} (\hat{u}_{xx} + \hat{u}_{yy}), \\ f_v &= \hat{v}_t + \lambda_{\text{adv}} (\hat{u}\,\hat{v}_x + \hat{v}\,\hat{v}_y) + \hat{p}_y - \lambda_{\text{vis}} (\hat{v}_{xx} + \hat{v}_{yy}), \\ f_c &= \hat{u}_x + \hat{v}_y , \end{align}\] and the residual loss \[\mathcal{L}_{r} \;=\; \frac{1}{N_f} \sum_{i=1}^{N_f} \Big( \,|f_u(t_i,x_i,y_i)|^2 \;+\; |f_v(t_i,x_i,y_i)|^2 \;+\; \alpha\,|f_c(t_i,x_i,y_i)|^2 \Big),\] where \(\alpha>0\) weights the continuity residual.

7.4.0.3 Data Loss (Vortex Shedding Velocities)

Given velocity measurements from the vortex-shedding flow, \(\{(t_i,x_i,y_i,u_i,v_i)\}_{i=1}^{N_d},\) we use only velocity supervision: \[\mathcal{L}_{d} \;=\; \frac{1}{N_d} \sum_{i=1}^{N_d} \Big( \,(\hat{u}(t_i,x_i,y_i) - u_i)^2 \;+\; (\hat{v}(t_i,x_i,y_i) - v_i)^2 \,\Big).\] Pressure is not supervised and is recovered implicitly via the momentum equations.

7.4.0.4 Total Loss

The training objective combines residual and data terms: \[\mathcal{L}(\theta,\lambda_{\text{adv}},\lambda_{\text{vis}}) \;=\; \lambda_{r}\,\mathcal{L}_{r} \;+\; \lambda_{d}\,\mathcal{L}_{d},\] with weights \(\lambda_r,\lambda_d>0\). For inverse problems, \(\lambda_{\text{adv}}\) and \(\lambda_{\text{vis}}\) are optimized jointly with \(\theta\); for forward problems they are fixed (e.g., \(\lambda_{\text{adv}}=1\), \(\lambda_{\text{vis}}=1/\mathrm{Re}\)).

8 Experimental Setup↩︎

8.1 Dataset and Task Construction↩︎

We evaluate our method across three PDE systems: 1D Convection, 2D Helmholtz, and 2D Navier-Stokes. For each PDE system, tasks are constructed by varying the governing PDE parameters following the experimental protocols established in prior PINN benchmark works.

8.1.0.1 1D Convection.

For the 1D Convection system, we use the analytical solution formulation introduced in Curr-Reg. [@krishnapriyan_characterizing_2021]. The dataset generation and task construction closely follow the publicly available benchmark implementation provided in their codebase. Task difficulty is controlled through the convection coefficient \(\beta\), where each task corresponds to a different value of \(\beta\).

8.1.0.2 2D Helmholtz.

For the 2D Helmholtz system, we construct our own dataset using the analytical solution formulation described in [@wang2021understanding]. The source term is analytically derived from the selected closed-form solution, allowing exact supervision of PDE behavior. Task construction is performed by varying the spatial frequency parameters \((a_1,a_2)\), which control the oscillatory complexity of the PDE solution.

8.1.0.3 2D Navier-Stokes.

For the 2D Navier-Stokes system, we generate fluid simulation data using PhiFlow, a differentiable fluid simulation framework based on finite element style numerical solvers. The simulation setup follows the cylinder wake configuration commonly adopted in PINN literature and is designed to closely match the benchmark setting introduced by [@raissi2017physicsinformeddeeplearning],. Task difficulty is controlled by varying the Reynolds number (\(Re\)), where higher Reynolds regimes correspond to increasingly complex flow dynamics.

Across all PDE systems, task families are constructed by systematically varying the corresponding PDE parameters, following the difficulty progression protocols established in the aforementioned baseline works.

8.2 Task Difficulty Characterization via PINN Loss Dynamics↩︎

In many PDE settings, easy vs. hard regimes are known and hence task difficulty is defined with the aid of domain expertise, where practitioners identify regimes that are known to be easy or challenging based on the underlying physics. However, in the absence of such prior knowledge, purely empirical criteria can be adopted, based on PINN training dynamics. Specifically, we observe that the distribution of PINN training losses even after early-stage training (i.e., 10K epochs) already provides a clear and consistent signal of task difficulty. Across PDE systems, this behavior exhibits a distinct transition in error dynamics, allowing us to separate regimes where standard PINN training is effective from those where it begins to fail.

8.2.0.1 1D Convection.

For the convection equation, task complexity is governed by the coefficient \(\beta\). As shown in Fig. 7 (a), lower values of \(\beta\) yield consistently low training loss across random seeds, indicating stable optimization. As \(\beta\) increases, the magnitude of the loss grows, with a sharp increase in magnitude at the transition from easy to difficult task boundary. Note, despite \(\beta=25\) exhibiting ‘hard task’ dynamics, it is included in the easy task set to maintain consistency with the experimental setup employed in Curr-Reg.( [@krishnapriyan_characterizing_2021]) a state-of-the-art baseline we compare with in our work. Based on this transition, we consider \(\beta \geq 30\) as the onset of the hard regime, and treat all tasks with \(\beta \geq 30\) as hard tasks in our experiments.

8.2.0.2 2D Helmholtz.

In the 2D Helmholtz equation, task complexity is controlled by the parameters \((a_1, a_2)\) which govern the domain spatial frequency (given by \(\sqrt{a_1^2 + a_2^2}\)). As shown in Fig. 7 (b), low-frequency configurations produce low and stable PINN training losses, indicating that PINNs can effectively learn task dynamics. In contrast, hard-task settings (i.e., high frequency configurations) exhibit a sharp increase in loss magnitude.

8.3 Baseline Model Description↩︎

We now provide brief descriptions of baseline models employed for comparative evaluation.

1. PINN (fixed) Raissi Original Formulation

2. PINN (dynamic) [@wang2021understanding; @wu2023comprehensive]: The current optimized and stable PINN training paradigm. In each iteration, collocation points are randomly resampled using quasi-random low-discrepancy sequences, combined with adaptive learning rate scheduling to improve convergence.

3. Curr-Reg. [@krishnapriyan_characterizing_2021]: Curriculum regularization of PINN which progressively trains models from easier to harder PDE tasks by incrementally introducing complexity during training. If we want to Learn a 1D convection system with \(\beta=30\). We incremently train on \(\beta=\{5,10,15,20,25,30\}\)

4. R3  [@daw_mitigating_2023]: Employs the adaptive sampling strategy Retain-Resample-Release (R3), focusing training effort by retaining collocation points with high residual errors while periodically resampling the remaining points.

5. Meta-PDE  [@qin2022metapdelearningsolvepdes]: A meta-learning framework for PDEs based on a MAML-style formulation. It also includes a learnable per-parameter step-size.

6. HyperLR [@cho2023hypernetwork]: A hypernetwork-based meta-learning method that generates low-rank, task-specific PINN weights for parameterized PDEs.

7. LIGO-PINN w/o (IE, GLO): Effectively standard PINN training from randomly initialized weights.

8. LIGO-PINN w/o GLO: An adaptation of the Model-Agnostic Meta-Learning–based framework.

8.4 Training Configuration per PDE System↩︎

We summarize the paper-facing experimental configurations used across all PDE systems. Unless otherwise specified, all models use Xavier weight initialization. Standard PINN-style baselines are optimized with Adam and cosine annealing learning-rate scheduling. For fair comparison, baseline methods and LIGO-PINN share the same network architectures, task definitions, collocation budgets, and evaluation protocols within each PDE system. The training horizon is matched to the method family. One-phase baselines are trained for 56,000 steps. LIGO/MAML-based downstream fine-tuning is run for 50,000 steps.

For LIGO-PINN, the LIGO invariance encoding stage is meta-trained exclusively on the easy-task distribution \(\Gamma_{\mathrm{Easy}}\), while evaluation is conducted on the unseen hard-task distribution \(\Gamma_{\mathrm{Hard}}\). During invariance encoding, the model learns shared initialization priors across related PDE tasks. It is then adapted to harder PDE regimes. Subsequently, gated layer-wise optimization (GLO) progressively adapts deeper network layers during task-specific fine-tuning. The corresponding outer and inner optimization steps and learning rates are reported explicitly in the tables below.

Unless explicitly stated otherwise, all reported MAE results are averaged across five random seeds. Within each PDE system, the same collocation point budgets, boundary-condition definitions, and evaluation protocols are used across competing methods. Table 5 reports the settings for the 1D Convection system, Table 6 reports the settings for the 2D Helmholtz system, and Table 7 reports the settings for the 2D Navier-Stokes system.

Reproducibility Notes. We provide the supplementary material needed to reproduce the reported results. Above, we specify the architecture, optimizer, learning-rate, task split, sampling, and training-horizon details used for each PDE system. The supplementary code package includes README-based instructions for environment setup, experiment execution, and result reproduction. We also provide the selected model runs and random-seed settings used for the reported evaluations so that the main benchmark tables and figures can be reproduced consistently.

Unless otherwise noted, all final reported MAE values are averaged over five random seeds. Analytical-task PDE systems are reproduced through the provided task-generation code, while the Navier–Stokes experiments use the supplied simulation/data-generation pipeline derived from the PhiFlow framework. Experiment metadata, metrics, and artifacts are tracked with MLflow for auditability and result verification.

8.5 Computational Resources↩︎

All experiments were conducted on an internal Linux GPU server using a four NVIDIA RTX A6000 GPU with 48 GB vRAM per GPU. To improve reproducibility, we report the hardware and software environment, the experiment execution setup, and representative runtime and memory costs for the main benchmark pipeline. The reported specifications correspond to the environment used for the paper-facing experiments.

8.5.0.1 Computation Cost.

To contextualize the computational overhead of invariance encoding, we report representative wall-clock time (minutes) and peak GPU memory usage (MB) in Table 9. All measurements were taken on the same NVIDIA RTX A6000 GPU to ensure consistency across stages. Wall-clock time was recorded per run, and peak GPU memory was measured using nvidia-smi. The table reports representative benchmark costs for the 1D Convection system and should be interpreted as stage-level execution cost rather than full cumulative project compute.

The overhead of LIGO-PINN is defined as the additional runtime introduced by the learned initialization meta-learning stage relative to PINN-D training. Across the three PDE systems, PINN-D requires 51.82 GPU-hours, while LIGO-PINN requires 62.57 GPU-hours when both learned initialization via meta-learning and downstream fine-tuning are included. This corresponds to an additional 10.75 GPU-hours, or a 20.75% cumulative runtime overhead. This overhead is incurred during the initialization stage and can be amortized when the learned initialization is reused across multiple downstream tasks within the same PDE family.

The higher peak memory observed for the Navier–Stokes system is primarily due to the larger dataset and the additional automatic-differentiation graph required by the coupled, higher-order PDE residual terms. Peak memory can be reduced by using a first-order MAML-style update, where higher-order gradients through the inner-loop adaptation are not tracked. This would improve runtime and memory efficiency, but it trades off meta-gradient fidelity and may reduce final accuracy. In this work, we keep the higher-order meta-gradient computation and use the same inner-loop setting across PDE systems to maintain experimental parity.

4pt

The invariance-encoding stage is a one-time cost. Once learned, the initialization can be reused across multiple tasks within the same PDE family, thereby amortizing the additional meta-training overhead. During the reported comparisons, the collocation budgets and evaluation settings were kept fixed between invariance encoding, downstream fine-tuning, and standard PINN training. This shows that LIGO-PINN retains comparable per-task fine-tuning cost to standard PINNs, with the main additional expense arising from the one-time initialization stage.

9 Additional Results↩︎

In this section, we present supplementary training results to further evaluate the behavior of our methods.

9.1 Do LIGO-PINN generalize effectively to PDEs defined on complex geometries?↩︎

To evaluate whether our method generalizes to PDEs on complex geometries, we compare LIGO-PINN with PINN (dynamic) for the screened Poisson equation defined on the Stanford Bunny surface [@Turk1994]. We further evaluate both PINN (dynamic) and LIGO-PINN under multiple random seed initializations to ensure robustness. Both models are trained for 50k epochs under identical training settings. LIGO-PINN achieves a 15% improvement over PINN (dynamic).

9.2 Hessian-Based Loss-Landscape Visualization↩︎

To better understand the local geometry of the trained PINN, we analyze the curvature of each loss component—total loss, data loss, and residual loss—around the converged weights \(\theta^*\). Concretely:

1. Calculate losses. Let \[L_{\mathrm{tot}}(\theta),\quad L_{\mathrm{data}}(\theta),\quad L_{\mathrm{res}}(\theta)\] denote the total, data, and residual losses, respectively, as in Eq. (2).

2. Compute Hessians. For each loss \(L\in\{L_{\mathrm{tot}},L_{\mathrm{data}},L_{\mathrm{res}}\}\), form the Hessian \[H_L = \nabla^2_\theta\,L(\theta)\bigl\lvert_{\theta=\theta^*}\,.\]

3. Eigen-decomposition. Solve \[H_L\,v_i = \lambda_i\,v_i,\quad \lambda_1 \ge \lambda_2 \ge \dots,\] and retain the top two eigenpairs \((\lambda_1,v_1)\), \((\lambda_2,v_2)\), which capture the directions of greatest curvature.

4. Parameter perturbation. For offsets \((\alpha,\beta)\in[-\delta,\delta]^2\), define \[\theta(\alpha,\beta) \;=\; \theta^* \;+\;\alpha\,v_1 \;+\;\beta\,v_2\,.\]

5. Landscape slice. The two-dimensional landscape \[f_L(\alpha,\beta) \;=\; L\bigl(\theta(\alpha,\beta)\bigr)\] is evaluated on a uniform grid of \((\alpha,\beta)\) and visualized as a heatmap.

2D Helmholtz Loss-Landscape Visualization. In 9 (a) we see that in case that model has converge to a trivial solution that has trapped in steep 1 dimension wells and stops which does nto faciliates learning. on the other hand 9 (b) we that a converged model has smoother optimization bowls that facilates good convergence

9.3 Exploring Effect of different initialization↩︎

Table 10 compares the extrapolation MAE for PINN-Dynamic under Xavier initialization, PINN with Kaiming initialization, and PINN-Dynamic with LIGO (i.e., nothing but our proposed LIGO-PINN model). The learned initialization (i.e., PINN-Dynamic LIGO) achieves consistently lower error than both Xavier and Kaiming baselines as task get progressively harder, indicating that the learned invariance weights provide a stronger starting point for optimization than standard random initializations.

4pt

Layer-wise Parameter Distribution Analysis. To characterize how the learned initialization modifies the PINN parameter landscape, we compare the layer-wise weight distributions of LIGO-PINN against layer-wise weight distributions of a randomly initialized PINN. We report (Table. 11) two statistics for the first, middle, and final layers across all PDEs namely: (i) the absolute difference between the variance of weight distributions per layer \(\Delta = |\mathrm{Var}(\text{LIGO-PINN})\)-\(\mathrm{Var}(\text{Random PINN})|\), and (ii) the 1-Wasserstein distance, which measures distributional similarity.

In Table 11, we find an interesting behavior where the results show a consistent pattern across systems. Specifically, we found a stark distributional dissimilarity (between randomly initialized PINNs and PINNs initialized with learned weights) in shallower layers with increased distributional similarity in deeper layers. This difference in distributional similarity can be attributed to a higher variance of weights in shallower layers of LIGO-PINN, compared to variance of weights in the corresponding layers of randomly initialized PINNs. Previous work has shown that neural networks with higher variance of weights in shallower layers have increased expressivity.[@yang2017mean], thereby allowing us to infer that LIGO-PINN with learned weights offer greater expressive power (and thereby larger learning capacity) compared to randomly initialized PINNs, also contributing to their superior performance in challenging PDE domains.

9.4 Compatibility with Existing PINN Stabilization Methods↩︎

We further investigate whether the proposed learned initialization strategy is compatible with existing PINN stabilization approaches. To this end, we integrate our learned initialization into the R3 sampling framework (R3) [@daw_mitigating_2023], a state-of-the-art adaptive collocation sampling method designed to improve PINN convergence stability.

Specifically, we compare two variants of the R3 framework:

• R3 (Xavier Init.): standard R3 training initialized using Xavier initialization,

• R3 (LIGO Init.): R3 initialized using learned weights obtained from our invariance encoding setup.

To provide a stronger stress test, both variants were trained for 100K epochs, significantly extending beyond the primary training budget used throughout this work. Evaluation was performed on challenging extrapolation regimes of the 1D Convection system with \(\beta \in \{40,50,60\}\). For fairness, both models used identical architectures, optimization schedules, and collocation settings. Results are averaged across 5 independent random seeds, and we report both mean MAE and standard deviation.

As shown in Table 12, the proposed learned initialization remains compatible with the R3 framework even under substantially extended training schedules. In particular, the learned initialization achieves lower MAE than Xavier initialization in the \(\beta=40\) and \(\beta=60\) extrapolation settings, including a substantial improvement in the more challenging \(\beta=60\) regime. These results suggest that the proposed initialization mechanism is complementary to existing PINN stabilization strategies and can be integrated into adaptive sampling frameworks in a model-agnostic manner.

9.5 Evaluation across baseline↩︎

9.5.0.1 1D Convection.

For the 1D convection system, we evaluate all methods under progressively harder extrapolation regimes with convection coefficients \(\beta \in \{40,50,60,70,80\}\). In addition to the mean prediction error, Table 13 also reports the variance across independent runs to characterize optimization stability under increasing transport difficulty. As the convection coefficient increases, standard PINN variants exhibit larger performance degradation and higher variance, indicating unstable convergence dynamics. In contrast, LIGO-PINN maintains consistently low error and reduced variance across most extrapolation settings, demonstrating improved robustness in convection-dominated regimes.

9.5.0.2 2D Helmholtz.

For the 2D Helmholtz system, we investigate extrapolation to high-frequency oscillatory regimes parameterized by \((a_1,a_2)\). These regimes introduce increasingly complex spatial oscillations that are known to challenge standard PINN optimization. Table 14 reports both the mean and variance of the MAE across independent runs to evaluate reconstruction quality and convergence stability. The results show that baseline PINN methods experience substantial instability and variance in higher-frequency settings, whereas LIGO-PINN consistently achieves lower error with significantly more stable convergence behavior across both isotropic and anisotropic wave propagation regimes.

9.5.0.3 2D Navier-Stokes (Forward Problem).

For the 2D incompressible Navier-Stokes system, we evaluate the forward pressure reconstruction problem across Reynolds numbers \(\mathrm{Re}\in\{600,800,1000\}\). In this setting, the pressure field is recovered without direct supervision, making optimization particularly challenging in high Reynolds number regimes with complex vortex dynamics. Table 15 reports both the mean and variance of the pressure MAE across independent runs to analyze convergence stability in this multi-physics setting. The results show that LIGO-PINN achieves substantially lower prediction error and reduced variability compared to existing baselines, indicating more stable and reliable convergence behavior in complex fluid dynamics regimes.

3pt

3pt

9.6 Ablation Analysis on 1D Convection and 2D Helmholtz↩︎

We extend the ablation analysis of LIGO-PINN to additional PDE domains, namely 1D Convection and 2D Helmholtz. As in the main text, we evaluate the effect of Invariance Encoding (IE) and Gated Layer-wise Optimization (GLO) using the ablation variants: ‘LIGO-PINN w/o (IE,GLO)’ (random initialization) and ‘LIGO-PINN w/o GLO’ (MAML-style initialization).

(i) Performance Trends. Across both PDE domains, random initialization (‘w/o (IE,GLO)’) results in higher error and larger variance, indicating unstable convergence behavior. Incorporating IE (‘w/o GLO’) improves performance and reduces variance relative to random initialization, demonstrating the benefit of task-aware initialization. However, the absence of GLO still leads to suboptimal convergence, particularly as task complexity increases. Full LIGO-PINN, combining IE and GLO, consistently achieves lower error and tighter error distributions, indicating improved training stability and robustness.

(ii) Spectral Bias Analysis. We analyze spectral bias by computing the Fourier transform of the absolute error fields. Across both 1D Convection and 2D Helmholtz, ablation variants exhibit higher spectral density, particularly in the high-frequency regime, indicating poor capture of fine-scale features. In contrast, LIGO-PINN consistently achieves lower spectral density across both low and high frequencies, demonstrating improved representation of complex solution structures. This confirms that the combination of invariance encoding and gated layer-wise optimization effectively mitigates spectral bias beyond what is achieved by either component alone.

9.7 Sensitivity Analysis of LIGO-PINN to Number and Range of IE Tasks↩︎

We further analyze the sensitivity of LIGO-PINN to two factors during the Invariance Encoding (IE) phase: (i) the number of sampled tasks \(k\), and (ii) the range of tasks used. For fairness, the baseline task ranges in the main paper were selected to align with prior work [@krishnapriyan_characterizing_2021], which we also compared against through the Curr-Reg. baseline.

9.7.0.1 Sensitivity to Number of Tasks \(k\).

We varied \(k \in \{2,3,4\}\) during meta-training, with \(k=3\) used in the main paper. Evaluation focused on the challenging extrapolation regimes of the 1D convection PDE, \(\beta \in \{50,60,70\}\). Table 16 reports the median MAE across random seeds. The results show that performance remains relatively stable across different values of \(k\), with \(k=3\) achieving the best median MAE on \(\beta=50\) and \(\beta=70\), while all settings remain in a similar range.

9.7.0.2 Sensitivity to Task Range.

We further test the robustness of LIGO-PINN by varying the range of easy tasks selected for IE. In addition to the main setup LIGO-PINN(5-25), we evaluate LIGO-PINN(10-30) and LIGO-PINN(15-35). The evaluation tasks remain the challenging extrapolation regimes \(\beta \in \{50,60,70\}\). Table 17 reports the median MAE across random seeds. The results show that all LIGO-PINN variants achieve comparable performance across task ranges, suggesting that the method is not overly sensitive to the exact IE task range.

9.7.0.3 Summary.

These results suggest that LIGO-PINN is robust to both (i) the number of sampled tasks \(k\) and (ii) the range of tasks selected during the IE phase. Across the tested settings, median MAE remains consistently low on challenging extrapolation tasks, indicating that the IE phase does not require highly precise tuning of either task count or task range. This supports the practical use of domain knowledge to select reasonable IE task ranges without sacrificing downstream extrapolation performance.