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

# Alternating Levenberg–Marquardt Training of Physics-Informed Neural Networks with Fourier-Enhanced Features
Source: [https://arxiv.org/html/2608.05892](https://arxiv.org/html/2608.05892)
###### Abstract

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 relativeL2L^\{2\}errors up to two orders of magnitude lower than state\-of\-the\-art baselines\.

###### keywords:

Physics\-informed neural network , Spectral bias , Levenberg–Marquardt algorithm , Alternating optimization

††journal:Journal of Computational Physics\\affiliation

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

## 1Introduction

Partial differential equations \(PDEs\) are fundamental tools for physical modeling, describing phenomena ranging from wave propagation and fluid flow to quantum dynamics\. Since analytical solutions are rarely available, classical methods are widely used to provide accurate approximations\. These include finite difference\[[1](https://arxiv.org/html/2608.05892#bib.bib1)\]and finite element methods\[[2](https://arxiv.org/html/2608.05892#bib.bib2),[3](https://arxiv.org/html/2608.05892#bib.bib3),[4](https://arxiv.org/html/2608.05892#bib.bib4)\], spectral methods\[[5](https://arxiv.org/html/2608.05892#bib.bib5),[6](https://arxiv.org/html/2608.05892#bib.bib6),[7](https://arxiv.org/html/2608.05892#bib.bib7)\], multi\-scale methods\[[8](https://arxiv.org/html/2608.05892#bib.bib8),[9](https://arxiv.org/html/2608.05892#bib.bib9),[10](https://arxiv.org/html/2608.05892#bib.bib10)\], and specialized approaches for highly oscillatory problems\[[11](https://arxiv.org/html/2608.05892#bib.bib11)\]\. Although accurate and well understood, these methods often rely on carefully constructed meshes or basis functions, which can be difficult to design for complex geometries, and their cost grows rapidly in high dimensions\. These limitations have motivated learning\-based solvers, including physics\-informed approaches that incorporate the governing equations directly into neural network training\.

Physics\-Informed Neural Networks \(PINNs\) have become one of the most widely used frameworks for solving PDEs\[[12](https://arxiv.org/html/2608.05892#bib.bib12),[13](https://arxiv.org/html/2608.05892#bib.bib13)\]\. A standard PINN represents the solution using a multi\-layer perceptron \(MLP\) and obtains its derivatives through automatic differentiation\. This enables the encoding of the governing equation, along with boundary and initial conditions, into a composite loss\. This mesh\-free formulation applies to both forward and inverse problems and extends naturally to challenging settings such as high\-dimensional PDEs\[[14](https://arxiv.org/html/2608.05892#bib.bib14)\], complex geometries\[[15](https://arxiv.org/html/2608.05892#bib.bib15)\], flow\-field reconstruction from sparse or noisy data\[[16](https://arxiv.org/html/2608.05892#bib.bib16),[17](https://arxiv.org/html/2608.05892#bib.bib17)\], and inverse design in nano\-optics and metamaterials\[[18](https://arxiv.org/html/2608.05892#bib.bib18)\]\. More broadly, physics\-informed learning has been integrated with Bayesian formulations for uncertainty quantification under noisy data\[[19](https://arxiv.org/html/2608.05892#bib.bib19)\], and with diffusion models that impose PDE constraints during generative sampling\[[20](https://arxiv.org/html/2608.05892#bib.bib20),[21](https://arxiv.org/html/2608.05892#bib.bib21),[22](https://arxiv.org/html/2608.05892#bib.bib22)\]\. It has also been extended to infinite\-dimensional, operator\-learning settings, such as Fourier Neural Operator\[[23](https://arxiv.org/html/2608.05892#bib.bib23)\]and DeepONet\[[24](https://arxiv.org/html/2608.05892#bib.bib24)\]\.

Despite this progress, several well\-known challenges still limit the effectiveness of PINNs\. A prominent issue isspectral bias, the tendency of neural networks to learn low\-frequency components first, which impedes the recovery of high\-frequency and multi\-scale solutions\[[25](https://arxiv.org/html/2608.05892#bib.bib25),[26](https://arxiv.org/html/2608.05892#bib.bib26)\]\. Training is also difficult for stiff, nonlinear, and chaotic problems, where standard PINNs may converge slowly or produce inaccurate solutions despite attaining a small training loss\[[27](https://arxiv.org/html/2608.05892#bib.bib27),[28](https://arxiv.org/html/2608.05892#bib.bib28)\]\. To address these challenges, a broad range of strategies has been developed to improve PINN trainability\.

Loss\-weighting methods adaptively balance the terms of the non\-convex composite objective\. For example, Wang et al\.\[[29](https://arxiv.org/html/2608.05892#bib.bib29)\]analyze the neural tangent kernel to equalize the convergence rates of the physics and data loss terms\. Pointwise weighting methods further adapt the contribution of individual points, loss\-attentional PINNs employ an auxiliary network that adaptively emphasizes high\-error points\[[30](https://arxiv.org/html/2608.05892#bib.bib30)\]; Anagnostopoulos et al\.\[[31](https://arxiv.org/html/2608.05892#bib.bib31)\]propose a gradient\-free, residual\-based attention scheme that updates pointwise weights directly from residual magnitudes; Si and Yan\[[32](https://arxiv.org/html/2608.05892#bib.bib32)\]extend pointwise weighting to local neighborhoods from a primal\-dual optimization perspective; and Zhao et al\.\[[33](https://arxiv.org/html/2608.05892#bib.bib33)\]impose temporal causality through non\-negative, monotonically decreasing pointwise weights that are decoupled from the collocation\-point arrangement\. A complementary direction improves the sampling of collocation points\. Wu et al\.\[[34](https://arxiv.org/html/2608.05892#bib.bib34)\]systematically compare non\-adaptive and residual\-based adaptive sampling strategies; Gao and Wang\[[35](https://arxiv.org/html/2608.05892#bib.bib35)\]select informative points by active learning for high\-dimensional nonlinear PDEs; and Lau et al\.\[[36](https://arxiv.org/html/2608.05892#bib.bib36)\]jointly select collocation and experimental points\.

Beyond reweighting and resampling, another line of work enriches the network representation itself through input embeddings, activation functions, and architectural modifications\. Fourier feature embeddings enrich the input representation with high\-frequency components\[[37](https://arxiv.org/html/2608.05892#bib.bib37)\], broadening the frequency spectrum accessible to the network therefore potentially alleviating spectral bias\[[38](https://arxiv.org/html/2608.05892#bib.bib38)\]\. At the level of activation functions, Sitzmann et al\.\[[39](https://arxiv.org/html/2608.05892#bib.bib39)\]adopt sinusoidal activations to represent oscillatory structures; Si et al\.\[[40](https://arxiv.org/html/2608.05892#bib.bib40)\]introduce learnable Cauchy activations to enhance the approximation of stiff and high\-frequency solutions; and Zeng and Zhu\[[41](https://arxiv.org/html/2608.05892#bib.bib41)\]introduce a sine\-enhanced adaptive activation for solving PDEs on irregular domains\. At the architectural level, Zhao et al\.\[[42](https://arxiv.org/html/2608.05892#bib.bib42)\]introduce transformer backbones that capture temporal dependencies in evolution equations; Wang et al\.\[[43](https://arxiv.org/html/2608.05892#bib.bib43)\]replace conventional MLP layers with Kolmogorov–Arnold layers built on learnable spline functions; Jagtap and Karniadakis\[[44](https://arxiv.org/html/2608.05892#bib.bib44)\]decompose the space\-time domain and assign a separate subnetwork to each subdomain; and Wang et al\.\[[45](https://arxiv.org/html/2608.05892#bib.bib45)\]stabilize the training of deep PINNs through adaptive residual connections combined with Fourier feature embeddings\.

However, accurately resolving highly oscillatory or strongly multi\-scale solutions remains challenging\. Fourier feature and activation\-based methods substantially alleviate spectral bias, but their accuracy may still plateau when the target solution spans a very wide frequency band or contains sharp structures, where even carefully tuned PINNs can struggle to achieve high precision\[[46](https://arxiv.org/html/2608.05892#bib.bib46),[47](https://arxiv.org/html/2608.05892#bib.bib47)\]\. More fundamentally, these advanced PINN variants train all network parameters jointly by minimizing a single nonconvex objective, in which feature learning and coefficient fitting remain entangled\. From an adaptive\-basis viewpoint, the effective basis changes throughout training while the output layer coefficients are simultaneously adjusted\[[48](https://arxiv.org/html/2608.05892#bib.bib48)\]\. This coupling can induce slow convergence, reduce robustness, and complicate theoretical analysis\.

Several methods mitigate this coupling by separating basis construction from coefficient fitting\. Extreme learning machines\[[49](https://arxiv.org/html/2608.05892#bib.bib49)\]and their physics\-informed variants\[[50](https://arxiv.org/html/2608.05892#bib.bib50),[51](https://arxiv.org/html/2608.05892#bib.bib51)\]fix a randomly sampled feature basis and recover the output weights by linear least squares, which makes the coefficient problem convex and efficient to solve\. However, the fixed basis limits approximation accuracy on complex or high\-frequency solutions\. Closest to our work, the iterative PINN with Fourier\-enhanced features \(IFeF\-PINN\)\[[52](https://arxiv.org/html/2608.05892#bib.bib52)\]makes the basis adaptive by applying a Fourier feature mapping to the last hidden layer, and is trained iteratively by updating the basis while solving for the output coefficients\. IFeF\-PINN is developed for scalar, linear PDEs, for which the coefficient subproblem is convex and admits a closed\-form solution; extending the decoupled framework to nonlinear PDE systems, where neither the closed\-form update nor the convexity\-based convergence analysis carries over, remains open\.

To bridge this gap, we propose the Fourier\-enhanced alternating Levenberg–Marquardt PINN \(FALM\-PINN\), a training framework that decouples adaptive basis learning from coefficient fitting and applies to general nonlinear and coupled PDE systems\. Training alternates between an upper\-level problem that updates the hidden layer parameters generating the Fourier\-enhanced basis, and a lower\-level problem that fits the projection coefficients on this basis using the Levenberg–Marquardt \(LM\) method\[[53](https://arxiv.org/html/2608.05892#bib.bib53),[54](https://arxiv.org/html/2608.05892#bib.bib54)\]\. Our contributions are threefold\. First, we formulate the lower\-level problem as a nonlinear least\-squares problem and solve it with the LM algorithm, which replaces the nonconvex problem with a sequence of strictly convex damped subproblems, each admitting a closed\-form solution\. The framework subsumes IFeF\-PINN\[[52](https://arxiv.org/html/2608.05892#bib.bib52)\]as a special case for linear PDEs\. Second, we show that the Fourier feature mapping applied to the learned latent representation induces an adaptive stationary kernel \(Lemma[1](https://arxiv.org/html/2608.05892#Thmlemma1)\), which interprets the adaptive basis learning as a form of deep kernel learning\[[55](https://arxiv.org/html/2608.05892#bib.bib55)\]and offers an explanation for its ability to mitigate spectral bias\. Third, we establish global convergence of the alternating scheme without requiring convexity of the lower\-level problem \(Theorem[1](https://arxiv.org/html/2608.05892#Thmtheorem1)\)\.

The organization of this paper is as follows\. Section[2](https://arxiv.org/html/2608.05892#S2)provides a brief overview of PINNs and Fourier feature mapping\. Section[3](https://arxiv.org/html/2608.05892#S3)introduces the FALM\-PINN framework, including the upper\- and lower\-level problems, the global alternating algorithm, and its convergence guarantees\. Section[4](https://arxiv.org/html/2608.05892#S4)reports numerical experiments on linear and nonlinear benchmarks, and Section[5](https://arxiv.org/html/2608.05892#S5)summarizes the conclusions and discusses some promising future directions\.

## 2Background

### 2\.1Physics\-informed neural networks

PINNs incorporate the governing equations directly into the training objective, enabling data\-driven learning under physical constraints\[[12](https://arxiv.org/html/2608.05892#bib.bib12),[13](https://arxiv.org/html/2608.05892#bib.bib13)\]\. Letn\>0n\>0denote the dimension of the spatio\-temporal domain,m\>0m\>0the number of solution components,M\>0M\>0the number of equations in the system, andk\>0k\>0the highest order of differentiation in the system\. We consider the following PDE system on a bounded domainΩ⊂ℝn\\Omega\\subset\\mathbb\{R\}^\{n\}:

𝒩​\[u\]​\(x\)\\displaystyle\\mathscr\{N\}\[u\]\(x\)=f​\(x\),x∈Ω,\\displaystyle=f\(x\),\\quad x\\in\\Omega,\(1\)ℬ​\[u\]​\(x\)\\displaystyle\\mathscr\{B\}\[u\]\(x\)=g​\(x\),x∈Γ⊆∂Ω,\\displaystyle=g\(x\),\\quad x\\in\\Gamma\\subseteq\\partial\\Omega,where𝒩:𝒰→C​\(Ω;ℝM\)\\mathscr\{N\}\\colon\\mathcal\{U\}\\to C\(\\Omega;\\mathbb\{R\}^\{M\}\)is a partial differential operator andℬ:𝒰→C​\(Γ;ℝm\)\\mathscr\{B\}\\colon\\mathcal\{U\}\\to C\(\\Gamma;\\mathbb\{R\}^\{m\}\)denotes the linear boundary operator, both defined on the classical solution space𝒰≔Ck​\(Ω¯;ℝm\)\\mathcal\{U\}\\coloneqq C^\{k\}\(\\overline\{\\Omega\};\\mathbb\{R\}^\{m\}\)\. This formulation accommodates a wide range of boundary conditions, including Dirichlet, Neumann, Robin, and periodic types\. The functionu∈𝒰u\\in\\mathcal\{U\}is the solution,f∈C​\(Ω;ℝM\)f\\in C\(\\Omega;\\mathbb\{R\}^\{M\}\)is the source term, andg∈C​\(Γ;ℝm\)g\\in C\(\\Gamma;\\mathbb\{R\}^\{m\}\)specifies the boundary condition\.

A standard PINN approximatesuuby an MLP, denoted byu^\\hat\{u\}, and is trained to minimize the composite loss:

ℒ^λ​\(u^\)=1Nb​∑i=1Nb‖ℬ​\[u^\]​\(xbi\)−g​\(xbi\)‖2\+λNf​∑i=1Nf‖𝒩​\[u^\]​\(xfi\)−f​\(xfi\)‖2,\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(\\hat\{u\}\)=\\frac\{1\}\{N\_\{b\}\}\\sum\_\{i=1\}^\{N\_\{b\}\}\\\|\\mathscr\{B\}\[\\hat\{u\}\]\(x\_\{b\}^\{i\}\)\-g\(x\_\{b\}^\{i\}\)\\\|^\{2\}\+\\frac\{\\lambda\}\{N\_\{f\}\}\\sum\_\{i=1\}^\{N\_\{f\}\}\\\|\\mathscr\{N\}\[\\hat\{u\}\]\(x\_\{f\}^\{i\}\)\-f\(x\_\{f\}^\{i\}\)\\\|^\{2\},\(2\)whereXb=\{xbi\}i=1Nb⊂ΓX\_\{b\}=\\\{x\_\{b\}^\{i\}\\\}\_\{i=1\}^\{N\_\{b\}\}\\subset\\GammaandXf=\{xfi\}i=1Nf⊂ΩX\_\{f\}=\\\{x\_\{f\}^\{i\}\\\}\_\{i=1\}^\{N\_\{f\}\}\\subset\\Omegaare the boundary and collocation point sets, respectively\. The physics weightλ\>0\\lambda\>0balances the two loss terms and may be fixed or adapted during training\. To explicitly illustrate the network architecture, consider anLL\-hidden layer network defined as:

z\(0\)\\displaystyle z^\{\(0\)\}=x,\\displaystyle=x,\(3\)z\(l\)\\displaystyle z^\{\(l\)\}=α​\(W\(l\)​z\(l−1\)\+c\(l\)\),1≤l≤L,\\displaystyle=\\alpha\\left\(W^\{\(l\)\}z^\{\(l\-1\)\}\+c^\{\(l\)\}\\right\),\\quad 1\\leq l\\leq L,u^​\(x\)\\displaystyle\\hat\{u\}\(x\)=W\(L\+1\)​z\(L\)\+c\(L\+1\),\\displaystyle=W^\{\(L\+1\)\}z^\{\(L\)\}\+c^\{\(L\+1\)\},whereα\\alphais the activation function \(e\.g\.,tanh\\tanh\),W\(l\)∈ℝdl×dl−1W^\{\(l\)\}\\in\\mathbb\{R\}^\{d\_\{l\}\\times d\_\{l\-1\}\}andc\(l\)∈ℝdlc^\{\(l\)\}\\in\\mathbb\{R\}^\{d\_\{l\}\}are the weight matrix and bias of layerll\. All hidden layer parameters are collected intoω=\{W\(l\),c\(l\)\}l=1L\\omega=\\\{W^\{\(l\)\},c^\{\(l\)\}\\\}\_\{l=1\}^\{L\}, and the final hidden layer output is denotedzω​\(x\)≜z\(L\)∈ℝpz\_\{\\omega\}\(x\)\\triangleq z^\{\(L\)\}\\in\\mathbb\{R\}^\{p\}, wherep=dLp=d\_\{L\}is the width of the final hidden layer\. This vector serves as a learned nonlinear latent representation of the input\. Finally,W\(L\+1\)∈ℝm×pW^\{\(L\+1\)\}\\in\\mathbb\{R\}^\{m\\times p\}andc\(L\+1\)∈ℝmc^\{\(L\+1\)\}\\in\\mathbb\{R\}^\{m\}are the output layer weight matrix and bias\. In standard PINN training, all parameters are optimized jointly by minimizingℒ^λ\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}via gradient descent\.

### 2\.2Fourier feature mapping

Random Fourier features, introduced by Rahimi and Recht\[[56](https://arxiv.org/html/2608.05892#bib.bib56)\], approximate stationary kernels via an explicit feature map and have since been adopted across machine learning and scientific computing\. The method rests on Bochner’s theorem, which characterizes any continuous, positive\-definite stationary kernel as the Fourier transform of a non\-negative measure\. Following Tancik et al\.\[[37](https://arxiv.org/html/2608.05892#bib.bib37)\], the Fourier feature mappingγD\\gamma\_\{D\}of an input coordinatev∈ℝqv\\in\\mathbb\{R\}^\{q\}is

γD​\(v\)=1D​\[cos⁡\(2​π​𝐁D​v\)sin⁡\(2​π​𝐁D​v\)\]∈ℝ2​D,\\gamma\_\{D\}\(v\)=\\frac\{1\}\{\\sqrt\{D\}\}\\begin\{bmatrix\}\\cos\(2\\pi\\mathbf\{B\}\_\{D\}v\)\\\\ \\sin\(2\\pi\\mathbf\{B\}\_\{D\}v\)\\end\{bmatrix\}\\in\\mathbb\{R\}^\{2D\},\(4\)where the entries of𝐁D∈ℝD×q\\mathbf\{B\}\_\{D\}\\in\\mathbb\{R\}^\{D\\times q\}are drawn independently from a Gaussian distribution𝒩​\(0,σ2\)\\mathcal\{N\}\(0,\\sigma^\{2\}\), and the bandwidthσ\>0\\sigma\>0controls the range of representable frequencies\. Training an MLP on these embedded coordinates corresponds to kernel regression with a stationary kernel of tunable bandwidth\[[37](https://arxiv.org/html/2608.05892#bib.bib37)\], which alleviates spectral bias and enables the network to learn high\-frequency components more efficiently\. The same mapping has been applied to the raw input coordinates in PINNs to improve the approximation of multi\-scale and high\-frequency PDEs\[[38](https://arxiv.org/html/2608.05892#bib.bib38)\]\.

## 3The FALM\-PINN framework

Standard PINN training updates the hidden layer parametersω\\omegaand the linear output layer jointly within a single nonconvex objective, entangling feature learning and coefficient fitting\. The latent featureszωz\_\{\\omega\}shift at every iteration while the output layer fits coefficients on a moving target\. This coupling, aggravated by the spectral bias that preventszωz\_\{\\omega\}from representing high\-frequency components, motivates training the two components separately\.

A natural idea is to decouple training into two alternating subproblems: an upper\-level problem that updates the hidden layer parametersω\\omegafor basis generation, and a lower\-level problem that solves for the projection coefficients\. This separation allows each problem to be addressed with a dedicated solver\. To further mitigate spectral bias, we apply the Fourier feature mapping of Section[2\.2](https://arxiv.org/html/2608.05892#S2.SS2)tozωz\_\{\\omega\}rather than to the raw input coordinates, enriching the basis with high\-frequency components on an adaptive feature space\. On this enriched basis, the lower\-level problem reduces to a regression that admits a closed\-form solution for linear PDEs and is solved efficiently by linearization for nonlinear ones\. We develop the two subproblems in Sections[3\.1](https://arxiv.org/html/2608.05892#S3.SS1)and[3\.2](https://arxiv.org/html/2608.05892#S3.SS2), combine them into the global alternating algorithm in Section[3\.3](https://arxiv.org/html/2608.05892#S3.SS3), and collect the theoretical properties of the framework in Section[3\.4](https://arxiv.org/html/2608.05892#S3.SS4)\. A fully worked example on a family of scalar conservation laws is provided in[A](https://arxiv.org/html/2608.05892#A1)\. The integrated architecture is illustrated in Fig\.[1](https://arxiv.org/html/2608.05892#S3.F1)\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/x1.png)Figure 1:Architecture of FALM\-PINN### 3\.1The upper\-level problem

The linear output layer in a standard PINN induces a dot\-product kernelk​\(x,x′\)=zω​\(x\)⊤​zω​\(x′\)k\(x,x^\{\\prime\}\)=z\_\{\\omega\}\(x\)^\{\\top\}z\_\{\\omega\}\(x^\{\\prime\}\)in the latent space, whose expressivity is determined entirely by the learned featureszωz\_\{\\omega\}\. To enrich this representation, we apply the Fourier feature mapping tozω​\(x\)z\_\{\\omega\}\(x\), replacing the dot\-product kernel with a stationary kernel that can represent high\-frequency components without increasing the network’s depth or width\.

Concretely, we apply the Fourier feature mappingγD:ℝp→ℝ2​D\\gamma\_\{D\}\\colon\\mathbb\{R\}^\{p\}\\to\\mathbb\{R\}^\{2D\}of Section[2\.2](https://arxiv.org/html/2608.05892#S2.SS2)to the latent representation, and define the composite feature mapψD:ℝn→ℝ2​D\\psi\_\{D\}\\colon\\mathbb\{R\}^\{n\}\\to\\mathbb\{R\}^\{2D\}as:

ψD​\(x\)=γD​\(zω​\(x\)\)=1D​\[cos⁡\(2​π​𝐁D​zω​\(x\)\)sin⁡\(2​π​𝐁D​zω​\(x\)\)\],\\psi\_\{D\}\(x\)=\\gamma\_\{D\}\\left\(z\_\{\\omega\}\(x\)\\right\)=\\frac\{1\}\{\\sqrt\{D\}\}\\begin\{bmatrix\}\\cos\\left\(2\\pi\\mathbf\{B\}\_\{D\}z\_\{\\omega\}\(x\)\\right\)\\\\ \\sin\\left\(2\\pi\\mathbf\{B\}\_\{D\}z\_\{\\omega\}\(x\)\\right\)\\end\{bmatrix\},\(5\)where𝐁D∈ℝD×p\\mathbf\{B\}\_\{D\}\\in\\mathbb\{R\}^\{D\\times p\}is a fixed Fourier feature matrix with entries drawn independently from𝒩​\(0,σ2\)\\mathcal\{N\}\(0,\\sigma^\{2\}\), andσ\>0\\sigma\>0is a bandwidth hyperparameter that governs the frequency range of the induced kernel\. We refer toψD​\(x\)\\psi\_\{D\}\(x\)as the Fourier\-enhanced features\. AsD→∞D\\to\\infty, the inner product⟨ψD​\(x\),ψD​\(x′\)⟩\\langle\\psi\_\{D\}\(x\),\\psi\_\{D\}\(x^\{\\prime\}\)\\rangleconverges almost surely to a Gaussian RBF kernel evaluated in the latent space \(Lemma[1](https://arxiv.org/html/2608.05892#Thmlemma1)\)\.

The approximated solution isuω,β=\[u1,…,um\]⊤u\_\{\\omega,\\beta\}=\\left\[u\_\{1\},\\dots,u\_\{m\}\\right\]^\{\\top\}withmmsolution components\. Projecting each component on the same Fourier\-enhanced basisψD​\(x\)∈ℝ2​D\\psi\_\{D\}\(x\)\\in\\mathbb\{R\}^\{2D\}, the final reconstructions will differ only in their projections’ coefficients\. Thekk\-th component is therefore parameterized asuk​\(x\)=⟨ψD​\(x\),θk⟩u\_\{k\}\(x\)=\\langle\\psi\_\{D\}\(x\),\\theta\_\{k\}\\ranglewith a trainable coefficient vectorθk∈ℝ2​D\\theta\_\{k\}\\in\\mathbb\{R\}^\{2D\}\. Collecting the coefficients column\-wise and stacking them into a single vector yields

Θ=\[θ1,…,θm\]∈ℝ2​D×m,\\displaystyle\\Theta=\[\\theta\_\{1\},\\ldots,\\theta\_\{m\}\]\\in\\mathbb\{R\}^\{2D\\times m\},\(6\)β=vec⁡\(Θ\)=\[θ1⊤,θ2⊤,…,θm⊤\]⊤∈ℝ2​m​D,\\displaystyle\\beta=\\operatorname\{vec\}\(\\Theta\)=\[\\theta\_\{1\}^\{\\top\},\\theta\_\{2\}^\{\\top\},\\ldots,\\theta\_\{m\}^\{\\top\}\]^\{\\top\}\\in\\mathbb\{R\}^\{2mD\},uω,β​\(x\)=Θ⊤​ψD​\(x\),x∈Ω\.\\displaystyle u\_\{\\omega,\\beta\}\(x\)=\\Theta^\{\\top\}\\psi\_\{D\}\(x\),\\quad x\\in\\Omega\.
Formally, with the coefficient vectorβ\\betaheld fixed, the upper\-level problem minimizes the shared training objective overω\\omegausing a gradient\-based optimizer \(e\.g\., Adam\[[57](https://arxiv.org/html/2608.05892#bib.bib57)\]\):

ω⋆​\(β\)=arg⁡minω⁡ℒ^λ​\(uω,β\)≔arg⁡minω⁡ℒupper​\(ω∣β\)\.\\omega^\{\\star\}\(\\beta\)=\\arg\\min\_\{\\omega\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\\bigl\(u\_\{\\omega,\\beta\}\\bigr\)\\coloneqq\\arg\\min\_\{\\omega\}\\mathcal\{L\}\_\{\\mathrm\{upper\}\}\(\\omega\\mid\\beta\)\.\(7\)The lower\-level problem minimizes the same objectiveℒ^λ\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}withω\\omegafixed, recast as an equivalent nonlinear least\-squares problem overβ\\beta\(Section[3\.2](https://arxiv.org/html/2608.05892#S3.SS2)\)\. The two updates are then alternated, as formalized in Section[3\.3](https://arxiv.org/html/2608.05892#S3.SS3)\.

### 3\.2The lower\-level problem

With the Fourier\-enhanced featuresψD\\psi\_\{D\}fixed for a givenω\\omega, the lower\-level problem finds the optimal projection coefficientsβ\\betathat satisfy the governing PDE and boundary conditions\. This decoupling exposes a least\-squares structure in the objective, which we exploit with the LM solver developed in Section[3\.2\.1](https://arxiv.org/html/2608.05892#S3.SS2.SSS1)\.

The PINN loss in Eq\. \([2](https://arxiv.org/html/2608.05892#S2.E2)\) has the form of a nonlinear least\-squares problem\. To set up the iterative solver of Section[3\.2\.1](https://arxiv.org/html/2608.05892#S3.SS2.SSS1), we define a composite residual vector\. Recall that the PDE system in Eq\. \([1](https://arxiv.org/html/2608.05892#S2.E1)\) comprisesMMequations andmmsolution components\. We construct the vectorized boundary data𝒚b∈ℝm​Nb\\bm\{y\}\_\{b\}\\in\\mathbb\{R\}^\{mN\_\{b\}\}and source term𝒇∈ℝM​Nf\\bm\{f\}\\in\\mathbb\{R\}^\{MN\_\{f\}\}by stacking over state variables and PDE equations, respectively:

𝒚b=\[𝒚b\(1\)⋮𝒚b\(m\)\],𝒇=\[𝒇\(1\)⋮𝒇\(M\)\],\\bm\{y\}\_\{b\}=\\begin\{bmatrix\}\\bm\{y\}\_\{b\}^\{\(1\)\}\\\\ \\vdots\\\\ \\bm\{y\}\_\{b\}^\{\(m\)\}\\end\{bmatrix\},\\quad\\bm\{f\}=\\begin\{bmatrix\}\\bm\{f\}^\{\(1\)\}\\\\ \\vdots\\\\ \\bm\{f\}^\{\(M\)\}\\end\{bmatrix\},\(8\)where each sub\-vector is defined as𝒚b\(k\)=\[gk​\(xb1\),…,gk​\(xbNb\)\]⊤∈ℝNb\\bm\{y\}\_\{b\}^\{\(k\)\}=\\big\[g\_\{k\}\(x\_\{b\}^\{1\}\),\\ldots,g\_\{k\}\(x\_\{b\}^\{N\_\{b\}\}\)\\big\]^\{\\top\}\\in\\mathbb\{R\}^\{N\_\{b\}\}fork=1,…,mk=1,\\dots,mand𝒇\(j\)=\[fj​\(xf1\),…,fj​\(xfNf\)\]⊤∈ℝNf\\bm\{f\}^\{\(j\)\}=\\big\[f\_\{j\}\(x\_\{f\}^\{1\}\),\\ldots,f\_\{j\}\(x\_\{f\}^\{N\_\{f\}\}\)\\big\]^\{\\top\}\\in\\mathbb\{R\}^\{N\_\{f\}\}forj=1,…,Mj=1,\\dots,M\.

Leveraging the linearity of the boundary operatorℬ\\mathscr\{B\}, we construct a matrix representation of the boundary residual𝑴b​\(ω\)∈ℝm​Nb×2​m​D\\bm\{M\}\_\{b\}\(\\omega\)\\in\\mathbb\{R\}^\{mN\_\{b\}\\times 2mD\}as:

𝑴b​\(ω\)=\[𝑴b\(1\)​\(ω\)⋮𝑴b\(m\)​\(ω\)\],𝑴b\(k\)​\(ω\)=\[ℬk​\[ψD\]​\(xb1\)⊤⋮ℬk​\[ψD\]​\(xbNb\)⊤\]∈ℝNb×2​m​D,\\bm\{M\}\_\{b\}\(\\omega\)=\\begin\{bmatrix\}\\bm\{M\}\_\{b\}^\{\(1\)\}\(\\omega\)\\\\ \\vdots\\\\ \\bm\{M\}\_\{b\}^\{\(m\)\}\(\\omega\)\\end\{bmatrix\},\\quad\\bm\{M\}\_\{b\}^\{\(k\)\}\(\\omega\)=\\begin\{bmatrix\}\\mathscr\{B\}\_\{k\}\[\\psi\_\{D\}\]\(x\_\{b\}^\{1\}\)^\{\\top\}\\\\ \\vdots\\\\ \\mathscr\{B\}\_\{k\}\[\\psi\_\{D\}\]\(x\_\{b\}^\{N\_\{b\}\}\)^\{\\top\}\\end\{bmatrix\}\\in\\mathbb\{R\}^\{N\_\{b\}\\times 2mD\},\(9\)whereℬk​\[ψD\]​\(xbi\)\\mathscr\{B\}\_\{k\}\[\\psi\_\{D\}\]\(x\_\{b\}^\{i\}\)represents the boundary operator applied to thekk\-th solution component, evaluated at pointxbix\_\{b\}^\{i\}, fori=1,…,Nbi=1,\\dots,N\_\{b\}\.

For the physics residual, the potential nonlinearity of the PDE operator𝒩\\mathscr\{N\}makes the mapping fromβ\\betato the residual nonlinear\. We collect the stacked physics terms into the nonlinear operator𝑴f​\(β∣ω\)∈ℝM​Nf\\bm\{M\}\_\{f\}\(\\beta\\mid\\omega\)\\in\\mathbb\{R\}^\{MN\_\{f\}\}, written in block\-vector form:

𝑴f​\(β∣ω\)=\[𝑴f\(1\)​\(β∣ω\)⋮𝑴f\(M\)​\(β∣ω\)\],𝑴f\(j\)​\(β∣ω\)=\[𝒩j​\[uω,β\]​\(xf1\)⋮𝒩j​\[uω,β\]​\(xfNf\)\]∈ℝNf,\\bm\{M\}\_\{f\}\(\\beta\\mid\\omega\)=\\begin\{bmatrix\}\\bm\{M\}\_\{f\}^\{\(1\)\}\(\\beta\\mid\\omega\)\\\\ \\vdots\\\\ \\bm\{M\}\_\{f\}^\{\(M\)\}\(\\beta\\mid\\omega\)\\end\{bmatrix\},\\quad\\bm\{M\}\_\{f\}^\{\(j\)\}\(\\beta\\mid\\omega\)=\\begin\{bmatrix\}\\mathscr\{N\}\_\{j\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{1\}\)\\\\ \\vdots\\\\ \\mathscr\{N\}\_\{j\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{N\_\{f\}\}\)\\end\{bmatrix\}\\in\\mathbb\{R\}^\{N\_\{f\}\},\(10\)where𝒩j​\[uω,β\]​\(xfi\)\\mathscr\{N\}\_\{j\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{i\}\)denotes thejj\-th equation of the PDE system evaluated atxfix\_\{f\}^\{i\}, fori=1,…,Nfi=1,\\dots,N\_\{f\}\. Substituting these formulations into the PINN objective yields the composite residual vector

𝒓​\(β∣ω\)=\[𝒓data​\(β∣ω\)𝒓phys​\(β∣ω\)\]∈ℝm​Nb\+M​Nf,\\bm\{r\}\(\\beta\\mid\\omega\)=\\begin\{bmatrix\}\\bm\{r\}\_\{\\mathrm\{data\}\}\(\\beta\\mid\\omega\)\\\\ \\bm\{r\}\_\{\\mathrm\{phys\}\}\(\\beta\\mid\\omega\)\\end\{bmatrix\}\\in\\mathbb\{R\}^\{mN\_\{b\}\+MN\_\{f\}\},\(11\)where

𝒓data​\(β∣ω\)\\displaystyle\\bm\{r\}\_\{\\mathrm\{data\}\}\(\\beta\\mid\\omega\)=1Nb​\(𝑴b​\(ω\)​β−𝒚b\)∈ℝm​Nb,\\displaystyle=\\frac\{1\}\{\\sqrt\{N\_\{b\}\}\}\\bigl\(\\bm\{M\}\_\{b\}\(\\omega\)\\beta\-\\bm\{y\}\_\{b\}\\bigr\)\\in\\mathbb\{R\}^\{mN\_\{b\}\},\(12\)𝒓phys​\(β∣ω\)\\displaystyle\\bm\{r\}\_\{\\mathrm\{phys\}\}\(\\beta\\mid\\omega\)=λNf​\(𝑴f​\(β∣ω\)−𝒇\)∈ℝM​Nf\.\\displaystyle=\\sqrt\{\\frac\{\\lambda\}\{N\_\{f\}\}\}\\bigl\(\\bm\{M\}\_\{f\}\(\\beta\\mid\\omega\)\-\\bm\{f\}\\bigr\)\\in\\mathbb\{R\}^\{MN\_\{f\}\}\.Under this setting, the lower\-level objective is a nonlinear least\-squares reformulation ofℒ^λ\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}overβ\\betafor any fixedω\\omega:

ℒlower\(β∣ω\)=∥𝒓\(β∣ω\)∥22\.\\mathcal\{L\}\_\{\\mathrm\{lower\}\}\(\\beta\\mid\\omega\)=\\\|\\bm\{r\}\(\\beta\\mid\\omega\)\\\|\_\{2\}^\{2\}\.\(13\)The optimal coefficient vectorβ⋆​\(ω\)\\beta^\{\\star\}\(\\omega\)is then obtained by minimizing:

β⋆​\(ω\)=arg⁡minβ⁡ℒ^λ​\(uω,β\)≔arg⁡minβ⁡ℒlower​\(β∣ω\)\.\\beta^\{\\star\}\(\\omega\)=\\arg\\min\_\{\\beta\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\\bigl\(u\_\{\\omega,\\beta\}\\bigr\)\\coloneqq\\arg\\min\_\{\\beta\}\\mathcal\{L\}\_\{\\mathrm\{lower\}\}\(\\beta\\mid\\omega\)\.\(14\)Since𝒓phys\\bm\{r\}\_\{\\mathrm\{phys\}\}is generally nonlinear inβ\\betathrough𝑴f​\(β∣ω\)\\bm\{M\}\_\{f\}\(\\beta\\mid\\omega\), the lower\-level objective is nonconvex and admits no closed\-form minimizer\. We propose to solve it iteratively in Section[3\.2\.1](https://arxiv.org/html/2608.05892#S3.SS2.SSS1)\.

#### 3\.2\.1Linearization via the Levenberg–Marquardt algorithm

We solve the lower\-level problem with the Levenberg–Marquardt \(LM\) algorithm\[[53](https://arxiv.org/html/2608.05892#bib.bib53),[54](https://arxiv.org/html/2608.05892#bib.bib54)\], which iteratively linearizes the residual and solves a sequence of damped least\-squares subproblems\. Since the upper\-level parametersω\\omegaremain fixed throughout the lower\-level solve, we write the residual vector and the stacked physics operator as𝒓​\(β\)≜𝒓​\(β∣ω\)\\bm\{r\}\(\\beta\)\\triangleq\\bm\{r\}\(\\beta\\mid\\omega\)and𝑴f​\(β\)≜𝑴f​\(β∣ω\)\\bm\{M\}\_\{f\}\(\\beta\)\\triangleq\\bm\{M\}\_\{f\}\(\\beta\\mid\\omega\)\.

##### First\-order linearization

Letβ\(j\)\\beta^\{\(j\)\}denote the current coefficient estimate at thejj\-th LM iteration\. We linearize the residual vector aroundβ\(j\)\\beta^\{\(j\)\}via a first\-order Taylor expansion:

𝒓​\(β\(j\)\+Δ​β\(j\)\)≈𝒓​\(β\(j\)\)\+𝐉j​Δ​β\(j\),\\displaystyle\\bm\{r\}\(\\beta^\{\(j\)\}\+\\Delta\\beta^\{\(j\)\}\)\\approx\\bm\{r\}\(\\beta^\{\(j\)\}\)\+\\mathbf\{J\}\_\{j\}\\Delta\\beta^\{\(j\)\},\(15\)whereΔ​β\(j\)∈ℝ2​m​D\\Delta\\beta^\{\(j\)\}\\in\\mathbb\{R\}^\{2mD\}is the update and𝐉j=∇β𝒓​\(β\)\|β=β\(j\)∈ℝ\(m​Nb\+M​Nf\)×2​m​D\\mathbf\{J\}\_\{j\}=\\nabla\_\{\\beta\}\\bm\{r\}\(\\beta\)\\big\|\_\{\\beta=\\beta^\{\(j\)\}\}\\in\\mathbb\{R\}^\{\(mN\_\{b\}\+MN\_\{f\}\)\\times 2mD\}is the Jacobian of the residual vector with the block form𝐉j=\[𝐉data⊤,𝐉phys\(j\)⊤\]⊤\\mathbf\{J\}\_\{j\}=\[\\mathbf\{J\}\_\{\\mathrm\{data\}\}^\{\\top\},\{\{\\mathbf\{J\}\_\{\\mathrm\{phys\}\}\}^\{\(j\)\}\}^\{\\top\}\]^\{\\top\},

𝐉data\\displaystyle\\mathbf\{J\}\_\{\\mathrm\{data\}\}=1Nb​𝑴b​\(ω\)∈ℝm​Nb×2​m​D,\\displaystyle=\\frac\{1\}\{\\sqrt\{N\_\{b\}\}\}\\bm\{M\}\_\{b\}\(\\omega\)\\in\\mathbb\{R\}^\{mN\_\{b\}\\times 2mD\},\(16\)𝐉phys\(j\)\\displaystyle\{\\mathbf\{J\}\_\{\\mathrm\{phys\}\}\}^\{\(j\)\}=λNf​∇β𝑴f​\(β\(j\)\)∈ℝM​Nf×2​m​D\.\\displaystyle=\\sqrt\{\\frac\{\\lambda\}\{N\_\{f\}\}\}\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)\\in\\mathbb\{R\}^\{MN\_\{f\}\\times 2mD\}\.Here,𝐉data\\mathbf\{J\}\_\{\\mathrm\{data\}\}is a constant matrix across iterations since𝒓data\\bm\{r\}\_\{\\mathrm\{data\}\}is linear inβ\\beta\.

##### Damped least\-squares subproblem

To iteratively minimize the nonlinear least\-squares objectiveℒlower​\(β∣ω\)\\mathcal\{L\}\_\{\\mathrm\{lower\}\}\(\\beta\\mid\\omega\)in Eq\. \([13](https://arxiv.org/html/2608.05892#S3.E13)\), we substitute the linearized residual \(LABEL:equ:lm\_taylor\) into the objective and augment with a damping term, obtaining the local quadratic subproblem:

minΔ​β\(j\)⁡Qj​\(Δ​β\(j\)\)=‖𝒓​\(β\(j\)\)\+𝐉j​Δ​β\(j\)‖22\+γj2​‖Δ​β\(j\)‖22,\\min\_\{\\Delta\\beta^\{\(j\)\}\}Q\_\{j\}\(\\Delta\\beta^\{\(j\)\}\)=\\bigl\\\|\\bm\{r\}\(\\beta^\{\(j\)\}\)\+\\mathbf\{J\}\_\{j\}\\Delta\\beta^\{\(j\)\}\\bigr\\\|\_\{2\}^\{2\}\+\\frac\{\\gamma\_\{j\}\}\{2\}\\\|\\Delta\\beta^\{\(j\)\}\\\|\_\{2\}^\{2\},\(17\)whereγj\>0\\gamma\_\{j\}\>0is the damping parameter at iterationjj\. Setting∇Δ​βQj=0\\nabla\_\{\\Delta\\beta\}Q\_\{j\}=0gives the closed\-form solution of the subproblem:

\(𝐉j⊤​𝐉j\+γj2​𝐈\)​Δ​β\(j\)=−𝐉j⊤​𝒓​\(β\(j\)\)\.\\bigl\(\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\+\\tfrac\{\\gamma\_\{j\}\}\{2\}\\mathbf\{I\}\\bigr\)\\Delta\\beta^\{\(j\)\}=\-\\mathbf\{J\}\_\{j\}^\{\\top\}\\bm\{r\}\(\\beta^\{\(j\)\}\)\.\(18\)The dynamic damping parameterγj\>0\\gamma\_\{j\}\>0ensures positive definiteness of the subproblem and balances the trade\-off between the Gauss–Newton and steepest descent directions\. For smallγj\\gamma\_\{j\}, the update approaches the Gauss–Newton step for fast convergence, whereΔ​β\(j\)≈−\(𝐉j⊤​𝐉j\)−1​𝐉j⊤​𝒓​\(β\(j\)\)\\Delta\\beta^\{\(j\)\}\\approx\-\(\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\)^\{\-1\}\\mathbf\{J\}\_\{j\}^\{\\top\}\\bm\{r\}\(\\beta^\{\(j\)\}\)\. In contrast, for largeγj\\gamma\_\{j\}, the step is dominated by the gradient, reducing to the steepest descent directionΔ​β\(j\)≈−\(2/γj\)​𝐉j⊤​𝒓​\(β\(j\)\)\\Delta\\beta^\{\(j\)\}\\approx\-\(2/\\gamma\_\{j\}\)\\mathbf\{J\}\_\{j\}^\{\\top\}\\bm\{r\}\(\\beta^\{\(j\)\}\), which guarantees a strict decrease inℒlower\\mathcal\{L\}\_\{\\mathrm\{lower\}\}\.

##### Summary of the lower\-level solver

Each lower\-level update solves a sequence ofJJdamped subproblems in Eq\. \([17](https://arxiv.org/html/2608.05892#S3.E17)\)\. The optimal projection coefficientsβ⋆​\(ω\)\\beta^\{\\star\}\(\\omega\), defined in Eq\. \([14](https://arxiv.org/html/2608.05892#S3.E14)\), are approximated by the LM outputβ\(J\)\\beta^\{\(J\)\}computed iteratively as:

Δ​β\(j\)\\displaystyle\\Delta\\beta^\{\(j\)\}=arg⁡minΔ​β⁡Qj​\(Δ​β\),\\displaystyle=\\arg\\min\_\{\\Delta\\beta\}Q\_\{j\}\(\\Delta\\beta\),\(19\)β\(j\+1\)\\displaystyle\\beta^\{\(j\+1\)\}=β\(j\)\+Δ​β\(j\),j=0,…,J−1,\\displaystyle=\\beta^\{\(j\)\}\+\\Delta\\beta^\{\(j\)\},\\quad j=0,\\dots,J\-1,with the initialβ\(0\)\\beta^\{\(0\)\}starting from the previous lower\-level solve\. Each subproblem is a strictly convex quadratic problem with a unique closed\-form solution given in Eq\. \([18](https://arxiv.org/html/2608.05892#S3.E18)\)\. Moreover, since the linearization error is quadratic in the residual \(Lemma[2](https://arxiv.org/html/2608.05892#Thmlemma2)\), the model is accurate with good initialization\. Empirically, one or two iterations \(J∈\{1,2\}J\\in\\\{1,2\\\}\) are sufficient to obtain an accurate estimate ofβ⋆​\(ω\)\\beta^\{\\star\}\(\\omega\)\. The per\-iteration cost is dominated by computing the Gram matrix𝐉j⊤​𝐉j\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}, which requires evaluating the physics Jacobian∇β𝑴f​\(β\(j\)\)\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)\. Since both the PDE operator𝒩\\mathscr\{N\}and the Fourier\-enhanced feature mapψD\\psi\_\{D\}are known in closed form, the Jacobian rows can be constructed analytically, as illustrated for the scalar conservation laws example in[A](https://arxiv.org/html/2608.05892#A1)\.

When the PDE operator𝒩\\mathscr\{N\}is linear inuu, the physics residual𝒓phys\\bm\{r\}\_\{\\mathrm\{phys\}\}becomes affine inβ\\beta, and therefore the linearization \(LABEL:equ:lm\_taylor\) is exact\. As the Jacobian𝐉\\mathbf\{J\}is independent ofβ\(j\)\\beta^\{\(j\)\}, the lower\-level problem reduces to a strictly convex Kernel Ridge Regression \(KRR\) problem whereℒlower\\mathcal\{L\}\_\{\\mathrm\{lower\}\}is a convex quadratic inβ\\beta\. Consequently, the LM iteration converges in a single step to its unique global minimizer\[[52](https://arxiv.org/html/2608.05892#bib.bib52)\]\.

### 3\.3Global alternating optimization

Combining the upper\- and lower\-level problems introduced above, we formulate the overall training as the following alternating optimization:

ω⋆​\(β\)\\displaystyle\\omega^\{\\star\}\(\\beta\)=arg⁡minω⁡ℒ^λ​\(uω,β\)≔arg⁡minω⁡ℒupper​\(ω∣β\),\\displaystyle=\\arg\\min\_\{\\omega\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\\bigl\(u\_\{\\omega,\\beta\}\\bigr\)\\coloneqq\\arg\\min\_\{\\omega\}\\mathcal\{L\}\_\{\\mathrm\{upper\}\}\(\\omega\\mid\\beta\),\(20\)β⋆​\(ω\)\\displaystyle\\beta^\{\\star\}\(\\omega\)=arg⁡minβ⁡ℒ^λ​\(uω,β\)≔arg⁡minβ⁡ℒlower​\(β∣ω\)\.\\displaystyle=\\arg\\min\_\{\\beta\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\\bigl\(u\_\{\\omega,\\beta\}\\bigr\)\\coloneqq\\arg\\min\_\{\\beta\}\\mathcal\{L\}\_\{\\mathrm\{lower\}\}\(\\beta\\mid\\omega\)\.A practical challenge is that the lower\-level solver requires both an expressive feature basisψD\\psi\_\{D\}and an initial coefficient estimate close toβ⋆​\(ω\)\\beta^\{\\star\}\(\\omega\)\. As shown in Lemma[2](https://arxiv.org/html/2608.05892#Thmlemma2), the LM linearization error scales with the squared residual, so a poor initialization yields inaccurate steps\. Whenω\\omegaandβ\\betaare both initialized randomly, the training may converge slowly or settle at a poor local minimum\. To address this, we precede the training with a joint warm\-up phase ofNwarmN\_\{\\mathrm\{warm\}\}iterations, during whichω\\omegaandβ\\betaare updated together on the objectiveℒ^λ\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}:

ωk\+1=ωk−ηω​∇ωℒ^λ​\(uωk,βk\),βk\+1=βk−ηβ​∇βℒ^λ​\(uωk,βk\)\.\\omega\_\{k\+1\}=\\omega\_\{k\}\-\\eta\_\{\\omega\}\\nabla\_\{\\omega\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(u\_\{\\omega\_\{k\},\\beta\_\{k\}\}\),\\qquad\\beta\_\{k\+1\}=\\beta\_\{k\}\-\\eta\_\{\\beta\}\\nabla\_\{\\beta\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(u\_\{\\omega\_\{k\},\\beta\_\{k\}\}\)\.\(21\)Algorithm[1](https://arxiv.org/html/2608.05892#alg1)summarizes the full training procedure, which consists of two phases\. Phase 1 performs the joint warm\-up of\(ω,β\)\(\\omega,\\beta\), providing an accurate initialization\. Phase 2 then applies the alternating optimization in Eq\. \([20](https://arxiv.org/html/2608.05892#S3.E20)\)\.

Algorithm 1FALM\-PINNInput:collocation sets

XbX\_\{b\},

XfX\_\{f\}; hyperparameters

λ\\lambda,

σ\\sigma,

DD,

JJ,

NwarmN\_\{\\mathrm\{warm\}\},

NalterN\_\{\\mathrm\{alter\}\},

ηω\\eta\_\{\\omega\},

ηβ\\eta\_\{\\beta\}
Initializeparameters

ω0\\omega\_\{0\},

β0\\beta\_\{0\}, and fixed Fourier feature matrix

𝐁D∼𝒩​\(0,σ2​I\)\\mathbf\{B\}\_\{D\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}I\)
⊳\\trianglerightPhase 1: Joint warm\-up

for

k=0k=0to

Nwarm−1N\_\{\\mathrm\{warm\}\}\-1do

Compute Fourier\-enhanced features

ψD​\(⋅\)=γD​\(zωk​\(⋅\)\)\\psi\_\{D\}\(\\cdot\)=\\gamma\_\{D\}\(z\_\{\\omega\_\{k\}\}\(\\cdot\)\)
Joint warm\-up:

ωk\+1=ωk−ηω​∇ωℒ^λ​\(uωk,βk\),βk\+1=βk−ηβ​∇βℒ^λ​\(uωk,βk\)\\omega\_\{k\+1\}=\\omega\_\{k\}\-\\eta\_\{\\omega\}\\nabla\_\{\\omega\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(u\_\{\\omega\_\{k\},\\beta\_\{k\}\}\),\\quad\\beta\_\{k\+1\}=\\beta\_\{k\}\-\\eta\_\{\\beta\}\\nabla\_\{\\beta\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(u\_\{\\omega\_\{k\},\\beta\_\{k\}\}\)
endfor

⊳\\trianglerightPhase 2: Alternating optimization

for

k=Nwarmk=N\_\{\\mathrm\{warm\}\}to

Nwarm\+Nalter−1N\_\{\\mathrm\{warm\}\}\+N\_\{\\mathrm\{alter\}\}\-1do

Compute Fourier\-enhanced features

ψD​\(⋅\)=γD​\(zωk​\(⋅\)\)\\psi\_\{D\}\(\\cdot\)=\\gamma\_\{D\}\(z\_\{\\omega\_\{k\}\}\(\\cdot\)\)
Lower\-level update:

Set inner solver initialization:

β\(0\)←βk\\beta^\{\(0\)\}\\leftarrow\\beta\_\{k\}
for

j=0j=0to

J−1J\-1do

Compute Jacobian

𝐉j\\mathbf\{J\}\_\{j\}via Eq\. \([16](https://arxiv.org/html/2608.05892#S3.E16)\)

Solve normal equation

\(𝐉j⊤​𝐉j\+γj2​𝐈\)​Δ​β\(j\)=−𝐉j⊤​𝒓​\(β\(j\)\)\\bigl\(\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\+\\tfrac\{\\gamma\_\{j\}\}\{2\}\\mathbf\{I\}\\bigr\)\\Delta\\beta^\{\(j\)\}=\-\\mathbf\{J\}\_\{j\}^\{\\top\}\\bm\{r\}\(\\beta^\{\(j\)\}\)
β\(j\+1\)=β\(j\)\+Δ​β\(j\)\\beta^\{\(j\+1\)\}=\\beta^\{\(j\)\}\+\\Delta\\beta^\{\(j\)\}

endfor

Set

βk\+1←β\(J\)\\beta\_\{k\+1\}\\leftarrow\\beta^\{\(J\)\}
Upper\-level update:

ωk\+1=ωk−ηω​∇ωℒ^λ​\(uωk,βk\+1\)\\omega\_\{k\+1\}=\\omega\_\{k\}\-\\eta\_\{\\omega\}\\nabla\_\{\\omega\}\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(u\_\{\\omega\_\{k\},\\beta\_\{k\+1\}\}\)
endfor

return

ωNwarm\+Nalter\\omega\_\{N\_\{\\mathrm\{warm\}\}\+N\_\{\\mathrm\{alter\}\}\},

βNwarm\+Nalter\\beta\_\{N\_\{\\mathrm\{warm\}\}\+N\_\{\\mathrm\{alter\}\}\}

### 3\.4Theoretical analysis

This section collects the theoretical properties of FALM\-PINN\. We first formalize the kernel interpretation of the Fourier\-enhanced basis, then bound the error of the LM linearization underlying the lower\-level solver, and finally establish global convergence of the alternating optimization\.

#### 3\.4\.1Kernel approximation in the latent space

Recall the composite feature mapψD=γD∘zω\\psi\_\{D\}=\\gamma\_\{D\}\\circ z\_\{\\omega\}of Section[3\.1](https://arxiv.org/html/2608.05892#S3.SS1)\. The following lemma formalizes its connection to a Gaussian kernel in the latent space\.

###### Lemma 1\(Kernel approximation\)\.

LetkRBF​\(z,z′\)=exp⁡\(−2​π2​σ2​‖z−z′‖22\)k\_\{\\mathrm\{RBF\}\}\(z,z^\{\\prime\}\)=\\exp\\left\(\-2\\pi^\{2\}\\sigma^\{2\}\\\|z\-z^\{\\prime\}\\\|\_\{2\}^\{2\}\\right\)be a Gaussian RBF kernel defined on the latent spaceℝp\\mathbb\{R\}^\{p\}\. For any spatio\-temporal inputsx,x′∈ℝnx,x^\{\\prime\}\\in\\mathbb\{R\}^\{n\}, letz=zω​\(x\)z=z\_\{\\omega\}\(x\)andz′=zω​\(x′\)z^\{\\prime\}=z\_\{\\omega\}\(x^\{\\prime\}\)\. Then the inner product of the explicit feature maps converges to the exact kernel almost surely:

limD→∞⟨ψD​\(x\),ψD​\(x′\)⟩=kRBF​\(z,z′\)\.\\lim\_\{D\\to\\infty\}\\langle\\psi\_\{D\}\(x\),\\psi\_\{D\}\(x^\{\\prime\}\)\\rangle=k\_\{\\mathrm\{RBF\}\}\(z,z^\{\\prime\}\)\.\(22\)

###### Proof\.

The inner product of the feature vectors is

⟨ψD​\(x\),ψD​\(x′\)⟩=1D​∑j=1D\[cos⁡\(2​π​𝐛j⊤​z\)​cos⁡\(2​π​𝐛j⊤​z′\)\+sin⁡\(2​π​𝐛j⊤​z\)​sin⁡\(2​π​𝐛j⊤​z′\)\],\\langle\\psi\_\{D\}\(x\),\\psi\_\{D\}\(x^\{\\prime\}\)\\rangle=\\frac\{1\}\{D\}\\sum\_\{j=1\}^\{D\}\\left\[\\cos\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}z\)\\cos\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}z^\{\\prime\}\)\+\\sin\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}z\)\\sin\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}z^\{\\prime\}\)\\right\],\(23\)where𝐛j⊤\\mathbf\{b\}\_\{j\}^\{\\top\}is thejj\-th row of𝐁D\\mathbf\{B\}\_\{D\}\. Applying the identitycos⁡α​cos⁡β\+sin⁡α​sin⁡β=cos⁡\(α−β\)\\cos\\alpha\\cos\\beta\+\\sin\\alpha\\sin\\beta=\\cos\(\\alpha\-\\beta\)yields:

⟨ψD​\(x\),ψD​\(x′\)⟩=1D​∑j=1Dcos⁡\(2​π​𝐛j⊤​\(z−z′\)\)\.\\langle\\psi\_\{D\}\(x\),\\psi\_\{D\}\(x^\{\\prime\}\)\\rangle=\\frac\{1\}\{D\}\\sum\_\{j=1\}^\{D\}\\cos\\left\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}\(z\-z^\{\\prime\}\)\\right\)\.\(24\)Since𝐛j∼𝒩​\(0,σ2​Ip\)\\mathbf\{b\}\_\{j\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}I\_\{p\}\), the expectation of each summand is the characteristic function of the Gaussian distribution evaluated at the displacementδ=z−z′\\delta=z\-z^\{\\prime\}:

𝔼​\[cos⁡\(2​π​𝐛j⊤​δ\)\]=Re⁡\(𝔼​\[ei​2​π​𝐛j⊤​δ\]\)=exp⁡\(−2​π2​σ2​‖δ‖22\)\.\\displaystyle\\mathbb\{E\}\\left\[\\cos\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}\\delta\)\\right\]=\\operatorname\{Re\}\\left\(\\mathbb\{E\}\\left\[e^\{i2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}\\delta\}\\right\]\\right\)=\\exp\\left\(\-2\\pi^\{2\}\\sigma^\{2\}\\\|\\delta\\\|\_\{2\}^\{2\}\\right\)\.\(25\)Since the𝐛j\\mathbf\{b\}\_\{j\}are i\.i\.d\. and\|cos⁡\(⋅\)\|≤1\|\\cos\(\\cdot\)\|\\leq 1, the strong law of large numbers applied to the sample mean in \([24](https://arxiv.org/html/2608.05892#S3.E24)\) gives:

1D​∑j=1Dcos⁡\(2​π​𝐛j⊤​δ\)→D→∞𝔼​\[cos⁡\(2​π​𝐛j⊤​δ\)\]=kRBF​\(z,z′\)almost surely\.\\frac\{1\}\{D\}\\sum\_\{j=1\}^\{D\}\\cos\\left\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}\\delta\\right\)\\xrightarrow\{D\\to\\infty\}\\mathbb\{E\}\\left\[\\cos\\left\(2\\pi\\mathbf\{b\}\_\{j\}^\{\\top\}\\delta\\right\)\\right\]=k\_\{\\mathrm\{RBF\}\}\(z,z^\{\\prime\}\)\\quad\\text\{almost surely\.\}∎

#### 3\.4\.2Accuracy of the LM linearization

We next quantify the accuracy of the first\-order linearization \(LABEL:equ:lm\_taylor\)\. The following bound shows that the linearization error is quadratic in the residual, which justifies the small number of inner iterations used in practice and motivates the warm\-up phase of Section[3\.3](https://arxiv.org/html/2608.05892#S3.SS3)\.

###### Lemma 2\(Linearization error\)\.

Assume the following:

1. 1\.𝒩\\mathscr\{N\}is twice continuously differentiable inuu\.
2. 2\.The lower\-level iterates lie in a compact setKKon which the residual Jacobian𝐉j\\mathbf\{J\}\_\{j\}has full column rank, with the smallest singular value bounded below byσ¯\>0\\underline\{\\sigma\}\>0\.
3. 3\.M𝒩<∞M\_\{\\mathscr\{N\}\}<\\inftyis a uniform bound overKKon the curvature of the PDE operator, with‖∇β2𝒩j​\[uω,β\]​\(xfi\)‖2≤M𝒩\\\|\\nabla^\{2\}\_\{\\beta\}\\mathscr\{N\}\_\{j\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{i\}\)\\\|\_\{2\}\\leq M\_\{\\mathscr\{N\}\}\.

Then the LM step \([18](https://arxiv.org/html/2608.05892#S3.E18)\) satisfies‖Δ​β\(j\)‖2≤σ¯−1​‖𝐫​\(β\(j\)\)‖2\\\|\\Delta\\beta^\{\(j\)\}\\\|\_\{2\}\\leq\\underline\{\\sigma\}^\{\-1\}\\\|\\bm\{r\}\(\\beta^\{\(j\)\}\)\\\|\_\{2\}, and the linearized residual \(LABEL:equ:lm\_taylor\) incurs the error

‖𝒓​\(β\(j\)\+Δ​β\(j\)\)−\[𝒓​\(β\(j\)\)\+𝐉j​Δ​β\(j\)\]‖2\\displaystyle\\bigl\\\|\\bm\{r\}\(\\beta^\{\(j\)\}\+\\Delta\\beta^\{\(j\)\}\)\-\\bigl\[\\bm\{r\}\(\\beta^\{\(j\)\}\)\+\\mathbf\{J\}\_\{j\}\\Delta\\beta^\{\(j\)\}\\bigr\]\\bigr\\\|\_\{2\}≤12​M𝒩​M​λ​‖Δ​β\(j\)‖22\\displaystyle\\leq\\tfrac\{1\}\{2\}M\_\{\\mathscr\{N\}\}\\sqrt\{M\\lambda\}\\,\\\|\\Delta\\beta^\{\(j\)\}\\\|\_\{2\}^\{2\}\(26\)≤M𝒩​M​λ2​σ¯2​‖𝒓​\(β\(j\)\)‖22\.\\displaystyle\\leq\\frac\{M\_\{\\mathscr\{N\}\}\\sqrt\{M\\lambda\}\}\{2\\underline\{\\sigma\}^\{2\}\}\\\|\\bm\{r\}\(\\beta^\{\(j\)\}\)\\\|\_\{2\}^\{2\}\.

###### Proof\.

The data residual is affine inβ\\beta, so its linearization is exact, and the error is therefore confined to theM​NfMN\_\{f\}physics components\. Each physics component has the formrphys,ℓ​\(β\)=λ/Nf​\(𝒩j​\[uω,β\]​\(xfi\)−fj​\(xfi\)\)r\_\{\\mathrm\{phys\},\\ell\}\(\\beta\)=\\sqrt\{\\lambda/N\_\{f\}\}\\,\(\\mathscr\{N\}\_\{j\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{i\}\)\-f\_\{j\}\(x\_\{f\}^\{i\}\)\)for some equation–point pair\(j,i\)\(j,i\)and isC2C^\{2\}inβ\\beta\. By the Taylor–Lagrange formula\[[58](https://arxiv.org/html/2608.05892#bib.bib58)\], there existsξℓ\\xi\_\{\\ell\}such that

\|rphys,ℓ​\(β\+Δ​β\)−rphys,ℓ​\(β\)−∇βrphys,ℓ​\(β\)⊤​Δ​β\|\\displaystyle\\bigl\|r\_\{\\mathrm\{phys\},\\ell\}\(\\beta\+\\Delta\\beta\)\-r\_\{\\mathrm\{phys\},\\ell\}\(\\beta\)\-\\nabla\_\{\\beta\}r\_\{\\mathrm\{phys\},\\ell\}\(\\beta\)^\{\\top\}\\Delta\\beta\\bigr\|=12​\|Δ​β⊤​∇β2rphys,ℓ​\(ξℓ\)​Δ​β\|\\displaystyle=\\tfrac\{1\}\{2\}\\bigl\|\\Delta\\beta^\{\\top\}\\nabla^\{2\}\_\{\\beta\}r\_\{\\mathrm\{phys\},\\ell\}\(\\xi\_\{\\ell\}\)\\,\\Delta\\beta\\bigr\|\(27\)≤12​λNf​M𝒩​‖Δ​β‖22\.\\displaystyle\\leq\\tfrac\{1\}\{2\}\\sqrt\{\\tfrac\{\\lambda\}\{N\_\{f\}\}\}\\,M\_\{\\mathscr\{N\}\}\\,\\\|\\Delta\\beta\\\|\_\{2\}^\{2\}\.Summing over theM​NfMN\_\{f\}physics components gives the first inequality in \([26](https://arxiv.org/html/2608.05892#S3.E26)\)\. The factor1/Nf1/\\sqrt\{N\_\{f\}\}in each component cancels theNf\\sqrt\{N\_\{f\}\}growth from theM​NfMN\_\{f\}components, leaving the bound independent ofNfN\_\{f\}\.

For the step bound, each subproblem \([17](https://arxiv.org/html/2608.05892#S3.E17)\) is a strictly convex quadratic and admits a closed\-form solution

Δ​β\(j\)=−\(𝐉j⊤​𝐉j\+γj2​𝐈\)−1​𝐉j⊤​𝒓​\(β\(j\)\)\.\\Delta\\beta^\{\(j\)\}=\-\\bigl\(\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\+\\tfrac\{\\gamma\_\{j\}\}\{2\}\\mathbf\{I\}\\bigr\)^\{\-1\}\\mathbf\{J\}\_\{j\}^\{\\top\}\\bm\{r\}\(\\beta^\{\(j\)\}\)\.\(28\)Writing the singular value decomposition \(SVD\)𝐉j=U​Σ​V⊤\\mathbf\{J\}\_\{j\}=U\\Sigma V^\{\\top\}with singular values\{σi\}\\\{\\sigma\_\{i\}\\\}, the operator in \([28](https://arxiv.org/html/2608.05892#S3.E28)\) can be written as

\(𝐉j⊤​𝐉j\+γj2​𝐈\)−1​𝐉j⊤=V​diag⁡\(σiσi2\+γj/2\)​U⊤,\\bigl\(\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\+\\tfrac\{\\gamma\_\{j\}\}\{2\}\\mathbf\{I\}\\bigr\)^\{\-1\}\\mathbf\{J\}\_\{j\}^\{\\top\}=V\\,\\operatorname\{diag\}\\\!\\Bigl\(\\tfrac\{\\sigma\_\{i\}\}\{\\sigma\_\{i\}^\{2\}\+\\gamma\_\{j\}/2\}\\Bigr\)U^\{\\top\},\(29\)whose singular values satisfy

σiσi2\+γj/2≤1σi≤1σ¯for any​γj≥0\.\\frac\{\\sigma\_\{i\}\}\{\\sigma\_\{i\}^\{2\}\+\\gamma\_\{j\}/2\}\\;\\leq\\;\\frac\{1\}\{\\sigma\_\{i\}\}\\;\\leq\\;\\frac\{1\}\{\\underline\{\\sigma\}\}\\qquad\\text\{for any \}\\gamma\_\{j\}\\geq 0\.\(30\)Taking the operator norm in \([28](https://arxiv.org/html/2608.05892#S3.E28)\) therefore gives

‖Δ​β\(j\)‖2≤σ¯−1​‖𝒓​\(β\(j\)\)‖2\.\\\|\\Delta\\beta^\{\(j\)\}\\\|\_\{2\}\\;\\leq\\;\\underline\{\\sigma\}^\{\-1\}\\,\\\|\\bm\{r\}\(\\beta^\{\(j\)\}\)\\\|\_\{2\}\.\(31\)Substituting into the first inequality of \([26](https://arxiv.org/html/2608.05892#S3.E26)\), evaluated atβ=β\(j\)\\beta=\\beta^\{\(j\)\}andΔ​β=Δ​β\(j\)\\Delta\\beta=\\Delta\\beta^\{\(j\)\}yields the second inequality\. Finally, if𝒩\\mathscr\{N\}is linear then∇β2𝒩j≡𝟎\\nabla^\{2\}\_\{\\beta\}\\mathscr\{N\}\_\{j\}\\equiv\\mathbf\{0\}, henceM𝒩=0M\_\{\\mathscr\{N\}\}=0and the error vanishes\. ∎

#### 3\.4\.3Convergence of the alternating optimization

We establish convergence guarantees for the alternating optimization in Phase 2 of Algorithm[1](https://arxiv.org/html/2608.05892#alg1)\. We denotez=\(ω,β\)z=\(\\omega,\\beta\)andΦ​\(z\):=ℒ^λ​\(uω,β\)\\Phi\(z\):=\\hat\{\\mathcal\{L\}\}\_\{\\lambda\}\(u\_\{\\omega,\\beta\}\)throughout this section, and analyze the sequence\{zk=\(ωk,βk\)\}\\\{z^\{k\}=\(\\omega^\{k\},\\beta^\{k\}\)\\\}generated by Phase 2\. We first state the key properties required for our convergence analysis\.

###### Assumption 1\(Bounded iterates and joint smoothness\)\.

The sequence\{zk\}k≥0\\\{z^\{k\}\\\}\_\{k\\geq 0\}remains in a compact set: there existsℬ\\mathcal\{B\}such thatzk∈ℬz^\{k\}\\in\\mathcal\{B\}for allk≥0k\\geq 0\.

SinceΦ\\Phiis a composition of smooth activation functions and a quadratic residual loss,Φ\\Phiis continuously differentiable onℬ\\mathcal\{B\}and∇2Φ\\nabla^\{2\}\\Phiis bounded there\. Consequently,Φ\\PhiisLL\-smooth onℬ\\mathcal\{B\}withL=supz∈ℬ‖∇2Φ​\(z\)‖2<∞L=\\sup\_\{z\\in\\mathcal\{B\}\}\\\|\\nabla^\{2\}\\Phi\(z\)\\\|\_\{2\}<\\infty, i\.e\.,

‖∇Φ​\(z\)−∇Φ​\(z′\)‖≤L​‖z−z′‖,∀z,z′∈ℬ\.\\\|\\nabla\\Phi\(z\)\-\\nabla\\Phi\(z^\{\\prime\}\)\\\|\\leq L\\\|z\-z^\{\\prime\}\\\|,\\qquad\\forall z,z^\{\\prime\}\\in\\mathcal\{B\}\.
###### Proposition 1\(Inexact LM descent condition\)\.

Every lower\-level update at iterationkkconsists ofJJdamped subproblems of the form Eq\. \([17](https://arxiv.org/html/2608.05892#S3.E17)\), withβ\(0\)=βk\\beta^\{\(0\)\}=\\beta^\{k\}andβ\(J\)=βk\+1\\beta^\{\(J\)\}=\\beta^\{k\+1\}\. There exist constantscβ,bβ\>0c\_\{\\beta\},b\_\{\\beta\}\>0such that the overall update satisfies:

Φ​\(ωk,βk\+1\)\\displaystyle\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)≤Φ​\(ωk,βk\)−cβ​‖βk\+1−βk‖2,\\displaystyle\\leq\\Phi\(\\omega^\{k\},\\beta^\{k\}\)\-c\_\{\\beta\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\},\(32\)‖∇βΦ​\(ωk,βk\+1\)‖\\displaystyle\\\|\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\\|≤bβ​‖βk\+1−βk‖\.\\displaystyle\\leq b\_\{\\beta\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|\.\(33\)

###### Proof\.

For brevity, we fixωk\\omega^\{k\}and writeΦ​\(⋅\)≡Φ​\(ωk,⋅\)\\Phi\(\\cdot\)\\equiv\\Phi\(\\omega^\{k\},\\cdot\)\. From the closed\-form LM update \([17](https://arxiv.org/html/2608.05892#S3.E17)\), the gradient identity holds:

∇βΦ​\(β\(j\)\)=2​𝐉j⊤​𝒓​\(β\(j\)\)=−\(2​𝐉j⊤​𝐉j\+γ​𝐈\)​Δ​β\(j\)\.\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(j\)\}\)=2\\mathbf\{J\}\_\{j\}^\{\\top\}\\bm\{r\}\(\\beta^\{\(j\)\}\)=\-\(2\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\+\\gamma\\mathbf\{I\}\)\\Delta\\beta^\{\(j\)\}\.
SinceΦ​\(β\)\\Phi\(\\beta\)isLL\-smooth from Assumption[1](https://arxiv.org/html/2608.05892#Thmassumption1),

Φ​\(β\(j\+1\)\)\\displaystyle\\Phi\(\\beta^\{\(j\+1\)\}\)≤Φ​\(β\(j\)\)\+⟨∇βΦ​\(β\(j\)\),Δ​β\(j\)⟩\+L2​‖Δ​β\(j\)‖2,\\displaystyle\\leq\\Phi\(\\beta^\{\(j\)\}\)\+\\langle\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(j\)\}\),\\Delta\\beta^\{\(j\)\}\\rangle\+\\tfrac\{L\}\{2\}\\\|\\Delta\\beta^\{\(j\)\}\\\|^\{2\},≤Φ​\(β\(j\)\)\+⟨−\(2​𝐉j⊤​𝐉j\+γ​𝐈\)​Δ​β\(j\),Δ​β\(j\)⟩\+L2​‖Δ​β\(j\)‖2,\\displaystyle\\leq\\Phi\(\\beta^\{\(j\)\}\)\+\\langle\-\(2\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}\+\\gamma\\mathbf\{I\}\)\\Delta\\beta^\{\(j\)\},\\Delta\\beta^\{\(j\)\}\\rangle\+\\tfrac\{L\}\{2\}\\\|\\Delta\\beta^\{\(j\)\}\\\|^\{2\},≤Φ​\(β\(j\)\)−2​‖𝐉j​Δ​β\(j\)‖2−\(γ−L2\)​‖Δ​β\(j\)‖2,\\displaystyle\\leq\\Phi\(\\beta^\{\(j\)\}\)\-2\\\|\\mathbf\{J\}\_\{j\}\\Delta\\beta^\{\(j\)\}\\\|^\{2\}\-\\big\(\\gamma\-\\tfrac\{L\}\{2\}\\big\)\\\|\\Delta\\beta^\{\(j\)\}\\\|^\{2\},≤Φ​\(β\(j\)\)−\(γ−L2\)​‖Δ​β\(j\)‖2\.\\displaystyle\\leq\\Phi\(\\beta^\{\(j\)\}\)\-\\big\(\\gamma\-\\tfrac\{L\}\{2\}\\big\)\\\|\\Delta\\beta^\{\(j\)\}\\\|^\{2\}\.With damping satisfyingγ\>L2\\gamma\>\\frac\{L\}\{2\}, we have the sufficient decrease condition for the single subproblem withc~β≔γ−L2\>0\\tilde\{c\}\_\{\\beta\}\\coloneqq\\gamma\-\\frac\{L\}\{2\}\>0\. Summing overj=0,1,…,J−1j=0,1,\\ldots,J\-1:

Φ​\(β\(J\)\)≤Φ​\(β\(0\)\)−c~β​∑j=0J−1‖β\(j\+1\)−β\(j\)‖2\.\\Phi\(\\beta^\{\(J\)\}\)\\leq\\Phi\(\\beta^\{\(0\)\}\)\-\\tilde\{c\}\_\{\\beta\}\\sum\_\{j=0\}^\{J\-1\}\\\|\\beta^\{\(j\+1\)\}\-\\beta^\{\(j\)\}\\\|^\{2\}\.By the Cauchy\-Schwarz inequality applied toβk\+1−βk=∑j=0J−1Δ​β\(j\)\\beta^\{k\+1\}\-\\beta^\{k\}=\\sum\_\{j=0\}^\{J\-1\}\\Delta\\beta^\{\(j\)\}:

‖βk\+1−βk‖2=‖∑j=0J−1Δ​β\(j\)‖2≤J​∑j=0J−1‖Δ​β\(j\)‖2,\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\}=\\big\\\|\\sum\_\{j=0\}^\{J\-1\}\\Delta\\beta^\{\(j\)\}\\big\\\|^\{2\}\\leq J\\sum\_\{j=0\}^\{J\-1\}\\\|\\Delta\\beta^\{\(j\)\}\\\|^\{2\},We have the sufficient decrease condition for the lower\-level problem as:

Φ​\(ωk,βk\+1\)≤Φ​\(ωk,βk\)−c~βJ​‖βk\+1−βk‖2,\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\leq\\Phi\(\\omega^\{k\},\\beta^\{k\}\)\-\\tfrac\{\\tilde\{c\}\_\{\\beta\}\}\{J\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\},Settingcβ:=c~β/Jc\_\{\\beta\}:=\\tilde\{c\}\_\{\\beta\}/Jproves Eq\. \([32](https://arxiv.org/html/2608.05892#S3.E32)\)\.

Next, we evaluate∇βΦ​\(β\(J\)\)\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\)\}\)by decomposition:

‖∇βΦ​\(β\(J\)\)‖\\displaystyle\\\|\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\)\}\)\\\|=‖∇βΦ​\(β\(J\)\)−∇βΦ​\(β\(J−1\)\)\+∇βΦ​\(β\(J−1\)\)‖\\displaystyle=\\\|\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\)\}\)\-\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\-1\)\}\)\+\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\-1\)\}\)\\\|≤‖∇βΦ​\(β\(J\)\)−∇βΦ​\(β\(J−1\)\)‖\+‖\(2​𝐉J−1⊤​𝐉J−1\+γ​𝐈\)​Δ​β\(J−1\)‖\\displaystyle\\leq\\\|\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\)\}\)\-\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\-1\)\}\)\\\|\+\\\|\(2\\mathbf\{J\}\_\{J\-1\}^\{\\top\}\\mathbf\{J\}\_\{J\-1\}\+\\gamma\\mathbf\{I\}\)\\Delta\\beta^\{\(J\-1\)\}\\\|≤L​‖Δ​β\(J−1\)‖\+\(2​‖𝐉J−1⊤​𝐉J−1‖2\+γ\)​‖Δ​β\(J−1\)‖,\\displaystyle\\leq L\\\|\\Delta\\beta^\{\(J\-1\)\}\\\|\+\\big\(2\\\|\\mathbf\{J\}\_\{J\-1\}^\{\\top\}\\mathbf\{J\}\_\{J\-1\}\\\|\_\{2\}\+\\gamma\\big\)\\\|\\Delta\\beta^\{\(J\-1\)\}\\\|,≤\(L\+2​‖𝐉J−1⊤​𝐉J−1‖2\+γ\)​‖Δ​β\(J−1\)‖\\displaystyle\\leq\\big\(L\+2\\\|\\mathbf\{J\}\_\{J\-1\}^\{\\top\}\\mathbf\{J\}\_\{J\-1\}\\\|\_\{2\}\+\\gamma\\big\)\\\|\\Delta\\beta^\{\(J\-1\)\}\\\|
Since the sequence\{zk\}k≥0\\\{z^\{k\}\\\}\_\{k\\geq 0\}is bounded on a compact set, we have‖𝐉J−1⊤​𝐉J−1‖2≤MJ\\\|\\mathbf\{J\}\_\{J\-1\}^\{\\top\}\\mathbf\{J\}\_\{J\-1\}\\\|\_\{2\}\\leq M\_\{J\}for some constantMJM\_\{J\}\. Definingb~β≔L\+2​MJ\+γ\\tilde\{b\}\_\{\\beta\}\\coloneqq L\+2M\_\{J\}\+\\gamma, we have‖∇βΦ​\(β\(J\)\)‖≤b~β​‖β\(J\)−β\(J−1\)‖\\\|\\nabla\_\{\\beta\}\\Phi\(\\beta^\{\(J\)\}\)\\\|\\leq\\tilde\{b\}\_\{\\beta\}\\\|\\beta^\{\(J\)\}\-\\beta^\{\(J\-1\)\}\\\|\. ForJ≥1J\\geq 1, there exists a finite constantκJ\>0\\kappa\_\{J\}\>0such that‖β\(J\)−β\(J−1\)‖≤κJ​‖∑j=0J−1Δ​β\(j\)‖=κJ​‖βk\+1−βk‖\\\|\\beta^\{\(J\)\}\-\\beta^\{\(J\-1\)\}\\\|\\leq\\kappa\_\{J\}\\big\\\|\\sum\_\{j=0\}^\{J\-1\}\\Delta\\beta^\{\(j\)\}\\big\\\|=\\kappa\_\{J\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|\. Settingbβ≔b~β​κJ\>0b\_\{\\beta\}\\coloneqq\\tilde\{b\}\_\{\\beta\}\\kappa\_\{J\}\>0proves Eq\. \([33](https://arxiv.org/html/2608.05892#S3.E33)\)\. ∎

###### Lemma 3\(Joint descent and gradient bound\)\.

With the same requirements as for Proposition[1](https://arxiv.org/html/2608.05892#Thmproposition1), if the upper\-level step size satisfiesηω<2/L\\eta\_\{\\omega\}<2/L, then there existsc=min⁡\(cβ,1ηω−L2\)\>0c=\\min\\left\(c\_\{\\beta\},\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\-\\tfrac\{L\}\{2\}\\right\)\>0andb=max⁡\(\(1ηω\+L\)2\+2​L2,bβ​2\)b=\\max\\left\(\\sqrt\{\\left\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\+L\\right\)^\{2\}\+2L^\{2\}\},b\_\{\\beta\}\\sqrt\{2\}\\right\)such that

Φ​\(zk\+1\)\\displaystyle\\Phi\(z^\{k\+1\}\)≤Φ​\(zk\)−c​‖zk\+1−zk‖2,\\displaystyle\\leq\\Phi\(z^\{k\}\)\-c\\\|z^\{k\+1\}\-z^\{k\}\\\|^\{2\},\(34\)‖∇Φ​\(zk\+1\)‖\\displaystyle\\\|\\nabla\\Phi\(z^\{k\+1\}\)\\\|≤b​‖zk\+1−zk‖\.\\displaystyle\\leq b\\\|z^\{k\+1\}\-z^\{k\}\\\|\.\(35\)

###### Proof\.

Each iterationkkconsists of two sub\-steps, withβ\\betaupdate\(ωk,βk\)→\(ωk,βk\+1\)\(\\omega^\{k\},\\beta^\{k\}\)\\to\(\\omega^\{k\},\\beta^\{k\+1\}\)through inner LM iterations, followed byω\\omegaupdate\(ωk,βk\+1\)→\(ωk\+1,βk\+1\)\(\\omega^\{k\},\\beta^\{k\+1\}\)\\to\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)through gradient descent\. Forβ\\betaupdate, Eq\. \([32](https://arxiv.org/html/2608.05892#S3.E32)\) directly gives

Φ​\(ωk,βk\+1\)≤Φ​\(ωk,βk\)−cβ​‖βk\+1−βk‖2\.\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\leq\\Phi\(\\omega^\{k\},\\beta^\{k\}\)\-c\_\{\\beta\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\}\.\(36\)
Forω\\omegaupdate, the descent lemma forLL\-smooth functions gives,

Φ​\(ωk\+1,βk\+1\)≤Φ​\(ωk,βk\+1\)\+⟨∇ωΦ​\(ωk,βk\+1\),ωk\+1−ωk⟩\+L2​‖ωk\+1−ωk‖2\.\\Phi\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)\\leq\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\+\\langle\\nabla\_\{\\omega\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\),\\omega^\{k\+1\}\-\\omega^\{k\}\\rangle\+\\tfrac\{L\}\{2\}\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\.Substituting the upper\-level update ruleωk\+1−ωk=−ηω​∇ωΦ​\(ωk,βk\+1\)\\omega^\{k\+1\}\-\\omega^\{k\}=\-\\eta\_\{\\omega\}\\nabla\_\{\\omega\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\), the inner\-product term becomes

⟨∇ωΦ​\(ωk,βk\+1\),ωk\+1−ωk⟩=−1ηω​‖ωk\+1−ωk‖2\.\\langle\\nabla\_\{\\omega\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\),\\omega^\{k\+1\}\-\\omega^\{k\}\\rangle=\-\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\.Hence we have

Φ​\(ωk\+1,βk\+1\)≤Φ​\(ωk,βk\+1\)−\(1ηω−L2\)​‖ωk\+1−ωk‖2\.\\Phi\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)\\leq\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\-\\big\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\-\\tfrac\{L\}\{2\}\\big\)\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\.\(37\)Definecω:=1ηω−L2c\_\{\\omega\}:=\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\-\\tfrac\{L\}\{2\}, for a small enough learning rateηω<2/L\\eta\_\{\\omega\}<2/Lwhich ensurescω\>0c\_\{\\omega\}\>0\. Combining Eq\. \([36](https://arxiv.org/html/2608.05892#S3.E36)\) and \([37](https://arxiv.org/html/2608.05892#S3.E37)\),

Φ​\(zk\+1\)≤Φ​\(zk\)−cβ​‖βk\+1−βk‖2−cω​‖ωk\+1−ωk‖2\.\\Phi\(z^\{k\+1\}\)\\leq\\Phi\(z^\{k\}\)\-c\_\{\\beta\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\}\-c\_\{\\omega\}\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\.Settingc:=min⁡\(cβ,cω\)\>0c:=\\min\(c\_\{\\beta\},c\_\{\\omega\}\)\>0proves Eq\. \([34](https://arxiv.org/html/2608.05892#S3.E34)\)\.

For the gradient bound, the update ofω\\omegaand jointLL\-smoothness of Assumption[1](https://arxiv.org/html/2608.05892#Thmassumption1)give

∇ωΦ​\(ωk,βk\+1\)=−1ηω​\(ωk\+1−ωk\),\\displaystyle\\nabla\_\{\\omega\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)=\-\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\(\\omega^\{k\+1\}\-\\omega^\{k\}\),\(38\)‖∇ωΦ​\(ωk\+1,βk\+1\)−∇ωΦ​\(ωk,βk\+1\)‖≤L​‖ωk\+1−ωk‖\.\\displaystyle\\\|\\nabla\_\{\\omega\}\\Phi\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)\-\\nabla\_\{\\omega\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\\|\\leq L\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|\.Combining both equations with the triangle inequality,

‖∇ωΦ​\(zk\+1\)‖≤\(1ηω\+L\)​‖ωk\+1−ωk‖\.\\\|\\nabla\_\{\\omega\}\\Phi\(z^\{k\+1\}\)\\\|\\leq\\big\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\+L\\big\)\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|\.\(39\)
For theβ\\betaupdate, we decompose the gradient as

‖∇βΦ​\(ωk\+1,βk\+1\)‖\\displaystyle\\\|\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)\\\|=‖\[∇βΦ​\(ωk\+1,βk\+1\)−∇βΦ​\(ωk,βk\+1\)\]\+∇βΦ​\(ωk,βk\+1\)‖,\\displaystyle=\\\|\\big\[\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)\-\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\big\]\+\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\\|,\(40\)≤‖∇βΦ​\(ωk\+1,βk\+1\)−∇βΦ​\(ωk,βk\+1\)‖\+‖∇βΦ​\(ωk,βk\+1\)‖,\\displaystyle\\leq\\\|\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\+1\},\\beta^\{k\+1\}\)\-\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\\|\+\\\|\\nabla\_\{\\beta\}\\Phi\(\\omega^\{k\},\\beta^\{k\+1\}\)\\\|,≤L​‖ωk\+1−ωk‖\+bβ​‖βk\+1−βk‖\.\\displaystyle\\leq L\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|\+b\_\{\\beta\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|\.Combining Eq\. \([39](https://arxiv.org/html/2608.05892#S3.E39)\) and \([40](https://arxiv.org/html/2608.05892#S3.E40)\),

‖∇Φ​\(zk\+1\)‖2\\displaystyle\\\|\\nabla\\Phi\(z^\{k\+1\}\)\\\|^\{2\}=‖∇ωΦ​\(zk\+1\)‖2\+‖∇βΦ​\(zk\+1\)‖2\\displaystyle=\\\|\\nabla\_\{\\omega\}\\Phi\(z^\{k\+1\}\)\\\|^\{2\}\+\\\|\\nabla\_\{\\beta\}\\Phi\(z^\{k\+1\}\)\\\|^\{2\}≤\(1ηω\+L\)2​‖ωk\+1−ωk‖2\+\(L​‖ωk\+1−ωk‖\+bβ​‖βk\+1−βk‖\)2,\\displaystyle\\leq\\big\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\+L\\big\)^\{2\}\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\+\\big\(L\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|\+\\ b\_\{\\beta\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|\\big\)^\{2\},≤\(1ηω\+L\)2​‖ωk\+1−ωk‖2\+2​L2​‖ωk\+1−ωk‖2\+2​bβ2​‖βk\+1−βk‖2\\displaystyle\\leq\\big\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\+L\\big\)^\{2\}\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\+2L^\{2\}\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\+2b\_\{\\beta\}^\{2\}\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\}=\(\(1ηω\+L\)2\+2​L2\)​\(‖ωk\+1−ωk‖2\)\+2​bβ2​\(‖βk\+1−βk‖2\)\.\\displaystyle=\\Big\(\\big\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\+L\\big\)^\{2\}\+2L^\{2\}\\Big\)\\big\(\\\|\\omega^\{k\+1\}\-\\omega^\{k\}\\\|^\{2\}\\big\)\+2b\_\{\\beta\}^\{2\}\\big\(\\\|\\beta^\{k\+1\}\-\\beta^\{k\}\\\|^\{2\}\\big\)\.Settingb=max⁡\(\(1ηω\+L\)2\+2​L2,bβ​2\)b=\\max\\left\(\\sqrt\{\\left\(\\tfrac\{1\}\{\\eta\_\{\\omega\}\}\+L\\right\)^\{2\}\+2L^\{2\}\},\\ b\_\{\\beta\}\\sqrt\{2\}\\right\)and taking square roots proves Eq\.\([35](https://arxiv.org/html/2608.05892#S3.E35)\)\.

∎

We now combine the joint sufficient decrease and gradient bound established in Lemma[3](https://arxiv.org/html/2608.05892#Thmlemma3)with the boundedness assumption to derive the convergence properties of the FALM\-PINN algorithm\.

###### Theorem 1\(Global convergence to critical points\)\.

The iterate sequence\{zk\}k≥0\\\{z^\{k\}\\\}\_\{k\\geq 0\}generated by Algorithm[1](https://arxiv.org/html/2608.05892#alg1)Phase 2 satisfies:

1. \(i\)\{Φ​\(zk\)\}k≥0\\\{\\Phi\(z^\{k\}\)\\\}\_\{k\\geq 0\}is non\-increasing and converges to some limitΦ⋆≥0\\Phi^\{\\star\}\\geq 0;
2. \(ii\)∑k=0∞‖zk\+1−zk‖2≤1c​\(Φ​\(z0\)−Φ⋆\)<∞\\sum\_\{k=0\}^\{\\infty\}\\\|z^\{k\+1\}\-z^\{k\}\\\|^\{2\}\\leq\\tfrac\{1\}\{c\}\\bigl\(\\Phi\(z^\{0\}\)\-\\Phi^\{\\star\}\\bigr\)<\\infty, and consequently‖zk\+1−zk‖→0\\\|z^\{k\+1\}\-z^\{k\}\\\|\\to 0ask→∞k\\to\\infty;
3. \(iii\)‖∇Φ​\(zk\)‖→0\\\|\\nabla\\Phi\(z^\{k\}\)\\\|\\to 0ask→∞k\\to\\infty;
4. \(iv\)\{zk\}\\\{z^\{k\}\\\}admits at least one accumulation point, and every accumulation pointz⋆z^\{\\star\}is a critical point ofΦ\\Phi, i\.e\.,∇Φ​\(z⋆\)=0\\nabla\\Phi\(z^\{\\star\}\)=0\.

###### Proof\.

The joint sufficient decrease in Lemma[3](https://arxiv.org/html/2608.05892#Thmlemma3)givesΦ​\(zk\+1\)≤Φ​\(zk\)−c​‖zk\+1−zk‖2≤Φ​\(zk\)\\Phi\(z^\{k\+1\}\)\\leq\\Phi\(z^\{k\}\)\-c\\\|z^\{k\+1\}\-z^\{k\}\\\|^\{2\}\\leq\\Phi\(z^\{k\}\)for everyk≥0k\\geq 0, so the sequence\{Φ​\(zk\)\}\\\{\\Phi\(z^\{k\}\)\\\}is non\-increasing\. SinceΦ​\(z\)≥0\\Phi\(z\)\\geq 0, the sequence is also bounded below by zero\. The monotone convergence theorem then yields a limitΦ⋆=limk→∞Φ​\(zk\)∈\[0,Φ​\(z0\)\]\\Phi^\{\\star\}=\\lim\_\{k\\to\\infty\}\\Phi\(z^\{k\}\)\\in\[0,\\Phi\(z^\{0\}\)\]\.

Rearranging the sufficient decrease inequality givesc​‖zk\+1−zk‖2≤Φ​\(zk\)−Φ​\(zk\+1\)c\\\|z^\{k\+1\}\-z^\{k\}\\\|^\{2\}\\leq\\Phi\(z^\{k\}\)\-\\Phi\(z^\{k\+1\}\), and summing overk=0,1,…,Nk=0,1,\\ldots,N:

c​∑k=0N‖zk\+1−zk‖2≤Φ​\(z0\)−Φ​\(zN\+1\)≤Φ​\(z0\)−Φ⋆\.c\\sum\_\{k=0\}^\{N\}\\\|z^\{k\+1\}\-z^\{k\}\\\|^\{2\}\\leq\\Phi\(z^\{0\}\)\-\\Phi\(z^\{N\+1\}\)\\leq\\Phi\(z^\{0\}\)\-\\Phi^\{\\star\}\.AsN→∞N\\to\\infty, the right\-hand side is a finite constant\. This implies that the infinite series is bounded, and hence we have‖zk\+1−zk‖→0\\\|z^\{k\+1\}\-z^\{k\}\\\|\\to 0\.

The joint gradient bound in Lemma[3](https://arxiv.org/html/2608.05892#Thmlemma3)gives‖∇Φ​\(zk\+1\)‖≤b​‖zk\+1−zk‖\\\|\\nabla\\Phi\(z^\{k\+1\}\)\\\|\\leq b\\\|z^\{k\+1\}\-z^\{k\}\\\|\. Then the right\-hand side converges to zero by \(ii\) ask→∞k\\to\\infty, thus‖∇Φ​\(zk\)‖→0\\\|\\nabla\\Phi\(z^\{k\}\)\\\|\\to 0\.

For the last claim, since the sequence\{zk\}⊂ℬ\\\{z^\{k\}\\\}\\subset\\mathcal\{B\}withℬ\\mathcal\{B\}compact from Assumption[1](https://arxiv.org/html/2608.05892#Thmassumption1), there exists a convergent subsequencezkj→z⋆∈ℬz^\{k\_\{j\}\}\\to z^\{\\star\}\\in\\mathcal\{B\}\. The jointLL\-smoothness in Assumption[1](https://arxiv.org/html/2608.05892#Thmassumption1)ensures that∇Φ\\nabla\\Phiis Lipschitz and therefore continuous, hence∇Φ​\(zkj\)→∇Φ​\(z⋆\)\\nabla\\Phi\(z^\{k\_\{j\}\}\)\\to\\nabla\\Phi\(z^\{\\star\}\)asj→∞j\\to\\infty\. By \(iii\),‖∇Φ​\(zkj\)‖→0\\\|\\nabla\\Phi\(z^\{k\_\{j\}\}\)\\\|\\to 0, and combining these two limits yields∇Φ​\(z⋆\)=0\\nabla\\Phi\(z^\{\\star\}\)=0\. ∎

## 4Numerical experiments

In this section, we conduct experiments to evaluate the proposed FALM\-PINN framework against representative and state\-of\-the\-art baselines across a range of PDE benchmarks, together with a plain PINN \(vanilla PINN\) as a reference\. Each baseline adopts a distinct strategy for mitigating spectral bias:

- 1\.Residual\-Based Attention \(RBA\)\[[31](https://arxiv.org/html/2608.05892#bib.bib31)\]: RBA is a gradient\-free, point\-wise reweighting strategy by calculating the moving average of normalized residuals, guiding the network to adaptively focus on high\-error regions\.
- 2\.Physics\-Informed Kolmogorov–Arnold Network \(PIKAN\)\[[43](https://arxiv.org/html/2608.05892#bib.bib43)\]: PIKAN leverages the Kolmogorov–Arnold Network architecture with learnable B\-spline functions to replace traditional MLPs, effectively mitigating spectral bias on high\-frequency and multi\-scale PDEs\.
- 3\.CompleX\-PINN\[[40](https://arxiv.org/html/2608.05892#bib.bib40)\]: CompleX\-PINN employs a learnable Cauchy activation function within a single hidden layer, reducing parameter complexity compared to conventional deep MLPs while targeting stiff and high\-frequency PDEs\.
- 4\.Sinusoidal Representation Networks \(SIREN\)\[[39](https://arxiv.org/html/2608.05892#bib.bib39)\]: SIREN replaces the standard smooth activation functions with periodic sinusoidal activationsϕi​\(x\)=sin⁡\(ω0​Wi​x\+bi\)\\phi\_\{i\}\(x\)=\\sin\(\\omega\_\{0\}W\_\{i\}x\+b\_\{i\}\), equipping the network with an inherent inductive bias toward high\-frequency functions\. This enables the network to represent fine\-scale solution features that are typically difficult to learn under spectral bias\.

##### Training and testing points

The testing dataset consists of uniformly distributed grid points: for 2D spatial or 1D spatio\-temporal PDEs, we use a uniform300×300300\\times 300grid of90,00090\{,\}000points, while for 2D spatio\-temporal PDEs, we use a uniform100×100×100100\\times 100\\times 100grid of100,000100\{,\}000points\. Training samples are drawn independently via Latin Hypercube Sampling \(LHS\) from the relevant domains, whereNfN\_\{f\},NbN\_\{b\}, andNiN\_\{i\}denote the number of physics collocation, boundary, and initial condition points, respectively\. The corresponding loss weights for these three terms are set to0\.010\.01,11, and11uniformly across all experiments\. The testing and training sets are not necessarily disjoint, as both are independently sampled from the same domain\.

Each method is trained and evaluated over five independent trials with different random seeds\. Unless specified, methods that employ a standard MLP use a fully connected network withtanh\\tanhactivations; methods with specialized architectures \(PIKAN, compleX\-PINN, SIREN\) follow their respective configurations\. All methods are optimized using Adam, and we follow the hyperparameter settings or adopt the tuning guidelines recommended for all baselines in their original publications\. Model accuracy is assessed using the relativeL2L^\{2\}error, defined as:

Relative​L2​error=∑k=1N\|u^​\(𝐱k\)−u​\(𝐱k\)\|2∑k=1N\|u​\(𝐱k\)\|2,\\text\{Relative \}L^\{2\}\\text\{ error\}=\\frac\{\\sqrt\{\\sum\\nolimits\_\{k=1\}^\{N\}\\bigl\|\\hat\{u\}\(\\mathbf\{x\}\_\{k\}\)\-u\(\\mathbf\{x\}\_\{k\}\)\\bigr\|^\{2\}\}\}\{\\sqrt\{\\sum\\nolimits\_\{k=1\}^\{N\}\\bigl\|u\(\\mathbf\{x\}\_\{k\}\)\\bigr\|^\{2\}\}\},\(41\)whereuudenotes the exact PDE solution,u^\\hat\{u\}is the predicted output of the model, andNNis the number of points in the test set\. We report the mean±\\pmstandard deviation across all five independent trials in the numerical tables, and present the heatmap results corresponding to the best\-performing trial for each method\. All experiments are conducted on a single NVIDIA RTX 4090 GPU\.

### 4\.12D Klein–Gordon equation

The Klein–Gordon equation serves as a representative test of long\-time approximation in nonlinear wave propagation\. We consider the following initial\-boundary value problem on the spatio\-temporal domainΩ=T×Ωs\\Omega=T\\times\\Omega\_\{s\}withT=\[0,10\]T=\[0,10\]andΩs=\[0,1\]2\\Omega\_\{s\}=\[0,1\]^\{2\}:

ut​t−Δ​u\+u2=f,\(t,x,y\)∈Ω,\\displaystyle u\_\{tt\}\-\\Delta u\+u^\{2\}=f,\\quad\(t,x,y\)\\in\\Omega,\(42\)u​\(0,x,y\)=x\+y,\(x,y\)∈Ωs,\\displaystyle u\(0,x,y\)=x\+y,\\quad\(x,y\)\\in\\Omega\_\{s\},ut​\(0,x,y\)=x​y,\(x,y\)∈Ωs,\\displaystyle u\_\{t\}\(0,x,y\)=xy,\\quad\(x,y\)\\in\\Omega\_\{s\},u​\(t,x,y\)=g​\(t,x,y\),\(x,y\)∈∂Ωs,t∈T,\\displaystyle u\(t,x,y\)=g\(t,x,y\),\\quad\(x,y\)\\in\\partial\\Omega\_\{s\},\\ t\\in T,where the source termffand boundary conditionggare obtained from the following analytical solution:

u​\(t,x,y\)=\(x\+y\)​cos⁡\(t\)\+x​y​sin⁡\(t\)\.u\(t,x,y\)=\(x\+y\)\\cos\(t\)\+xy\\sin\(t\)\.\(43\)
All MLP\-based methods use a fully connected network with44hidden layers of width8080, and are trained for50,00050\{,\}000iterations withNf=1,000N\_\{f\}=1\{,\}000collocation points,Nb=300N\_\{b\}=300boundary points, andNi=100N\_\{i\}=100initial points\. The baselines use Adam with an initial learning rate5×10−35\\times 10^\{\-3\}, decayed by a factor of0\.80\.8every1,0001\{,\}000iterations to a minimum of1×10−51\\times 10^\{\-5\}\.

Following Algorithm[1](https://arxiv.org/html/2608.05892#alg1), training is split into two phases\. In the warm\-up phase \(Nwarm=5,000N\_\{\\mathrm\{warm\}\}=5\{,\}000iterations\), bothω\\omegaandβ\\betaare jointly optimized by Adam with an initial learning rate of5×10−35\\times 10^\{\-3\}, decayed by a factor of0\.80\.8every1,0001\{,\}000iterations\. The Fourier feature mapping usesD=800D=800features with bandwidthσ=1\\sigma=1\. In the subsequent alternating optimization phase \(Nalter=45,000N\_\{\\mathrm\{alter\}\}=45\{,\}000iterations\), the upper\-level parametersω\\omegaare updated by Adam with a reduced learning rate of1×10−41\\times 10^\{\-4\}, following the same decay schedule to a floor of1×10−51\\times 10^\{\-5\}\. Each lower\-level solve performsJ=2J=2damped subproblems with a fixed damping parameterγ=1×10−6\\gamma=1\\times 10^\{\-6\}\.

The method\-specific hyperparameters for each baseline are as follows\. The RBA approach\[[31](https://arxiv.org/html/2608.05892#bib.bib31)\]applies residual\-based attention reweighting with parametersηRBA=0\.001\\eta\_\{\\mathrm\{RBA\}\}=0\.001andγRBA=0\.999\\gamma\_\{\\mathrm\{RBA\}\}=0\.999\. CompleX\-PINN\[[40](https://arxiv.org/html/2608.05892#bib.bib40)\]employs a single Cauchy layer of width200200\. PIKAN\[[43](https://arxiv.org/html/2608.05892#bib.bib43)\]uses a Kolmogorov–Arnold Network of hidden structure\[3,8,8,8,1\]\[3,8,8,8,1\], grid size1010, and spline order33\. SIREN\[[39](https://arxiv.org/html/2608.05892#bib.bib39)\]employs sinusoidal activations following the original implementation, withω0=30\\omega\_\{0\}=30\.

As illustrated in Fig\.[2](https://arxiv.org/html/2608.05892#S4.F2), FALM\-PINN exhibits a distinct two\-stage behavior: during the initial warm\-up phase, its convergence rate remains comparable to that of the baselines\. After switching to the alternating optimization phase, the error decreases by nearly two orders of magnitude while maintaining stable convergence throughout training\. This sharp improvement indicates that the warm\-up stage provides a suitable feature basis and coefficient initialization, while the lower\-level subproblem updates in Eq\. \([17](https://arxiv.org/html/2608.05892#S3.E17)\) rapidly refine the projection coefficients to produce a more accurate approximation\. This demonstrates the advantage of decoupling basis learning from coefficient optimization for nonlinear PDEs\. The final average relativeL2L^\{2\}errors, reported in Table[1](https://arxiv.org/html/2608.05892#S4.T1), show that FALM\-PINN achieves errors two orders of magnitude lower than all baselines\. Meanwhile, its training time and peak memory remain within the range spanned by the baselines, indicating that this substantial accuracy gain is not achieved at the expense of excessive computational cost\.

Table 1:2D Klein–Gordon equation: relativeL2L^\{2\}errors \(mean±\\pmstd over five independent runs\), GPU time per10001000iterations, and peak GPU memory for all methods\. The lowest value in each column is shown in bold\.MethodRelativeL2L^\{2\}errorTime \(s\)Memory \(GB\)Vanilla PINN8\.11×10−3±3\.98×10−48\.11\\times 10^\{\-3\}\\pm 3\.98\\times 10^\{\-4\}4\.06\\mathbf\{4\.06\}0\.07\\mathbf\{0\.07\}RBA2\.96×10−3±1\.05×10−32\.96\\times 10^\{\-3\}\\pm 1\.05\\times 10^\{\-3\}4\.074\.070\.640\.64compleX\-PINN4\.03×10−3±4\.94×10−44\.03\\times 10^\{\-3\}\\pm 4\.94\\times 10^\{\-4\}5\.055\.053\.773\.77PIKAN1\.14×10−3±1\.68×10−41\.14\\times 10^\{\-3\}\\pm 1\.68\\times 10^\{\-4\}36\.0336\.031\.951\.95SIREN4\.02×10−3±7\.79×10−44\.02\\times 10^\{\-3\}\\pm 7\.79\\times 10^\{\-4\}5\.815\.810\.950\.95FALM\-PINN3\.50×𝟏𝟎−𝟓±4\.88×𝟏𝟎−𝟔\\mathbf\{3\.50\\times 10^\{\-5\}\\pm 4\.88\\times 10^\{\-6\}\}7\.947\.940\.250\.25![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/comparison_kg2d_plot.png)Figure 2:Convergence trajectories on the 2D Klein–Gordon equation\. RelativeL2L^\{2\}error against training iterations; curves denote the mean across five independent trials, and shaded bands indicate the min–max range\. The vertical dashed line marks the transition from the warm\-up phase \(Nwarm=5,000N\_\{\\mathrm\{warm\}\}=5\{,\}000\) to the alternating optimization phase\.To further assess the approximation, Fig\.[3](https://arxiv.org/html/2608.05892#S4.F3)and[4](https://arxiv.org/html/2608.05892#S4.F4)display the exact solution and point\-wise absolute error maps for all methods at three time slicest∈\{2\.5,5\.0,7\.5\}t\\in\\\{2\.5,5\.0,7\.5\\\}\. FALM\-PINN yields consistently low errors throughout the spatial domain and over long\-time evolution, with the maximum point\-wise absolute error remaining below3×10−53\\times 10^\{\-5\}\. These results validate the ability of FALM\-PINN to accurately resolve nonlinear PDEs over long time horizons with a small training set\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/kg2d.png)Figure 3:Exact solution of the 2D Klein–Gordon equation at three time slicest∈\{2\.5,5\.0,7\.5\}t\\in\\\{2\.5,5\.0,7\.5\\\}\.![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/kg2d_solution_comparison.png)Figure 4:Point\-wise absolute error maps in log scale on the 2D Klein–Gordon equation at time slicest∈\{2\.5,5\.0,7\.5\}t\\in\\\{2\.5,5\.0,7\.5\\\}\(columns\) for vanilla PINN, RBA, compleX\-PINN, PIKAN, SIREN, and FALM\-PINN respectively\.
### 4\.21D Korteweg–de Vries equation

The Korteweg–de Vries \(KdV\) equation couples nonlinear convection with third\-order dispersion, producing dispersive wave dynamics that pose a challenge for PINNs\. Following\[[59](https://arxiv.org/html/2608.05892#bib.bib59)\], we consider the equation over the short time horizon on the spatio\-temporal domainΩ=T×Ωs\\Omega=T\\times\\Omega\_\{s\}withΩs=\[−1,1\]\\Omega\_\{s\}=\[\-1,1\]andT=\[0,1\]T=\[0,1\]\.

ut\+u​ux\+0\.0025​ux​x​x=0,\(t,x\)∈Ω,\\displaystyle u\_\{t\}\+uu\_\{x\}\+0025\\,u\_\{xxx\}=0,\\quad\(t,x\)\\in\\Omega,\(44\)u​\(0,x\)=cos⁡\(π​x\),x∈Ωs,\\displaystyle u\(0,x\)=\\cos\(\\pi x\),\\quad x\\in\\Omega\_\{s\},u​\(t,−1\)=u​\(t,1\),ux​\(t,−1\)=ux​\(t,1\),ux​x​\(t,−1\)=ux​x​\(t,1\),t∈T\.\\displaystyle u\(t,\-1\)=u\(t,1\),\\ \\ u\_\{x\}\(t,\-1\)=u\_\{x\}\(t,1\),\\ \\ u\_\{xx\}\(t,\-1\)=u\_\{xx\}\(t,1\),\\quad t\\in T\.The reference solution is shown in Fig\.[5](https://arxiv.org/html/2608.05892#S4.F5)\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/kdv_exact.png)Figure 5:Reference solution of the 1D Korteweg–de Vries equation\.All MLP\-based methods share the same architecture of66hidden layers of width6464\. We useNf=10,000N\_\{f\}=10\{,\}000collocation points,Ni=Nb=200N\_\{i\}=N\_\{b\}=200initial and boundary points\. The baselines are optimized by Adam with an initial learning rate of2×10−32\\times 10^\{\-3\}, exponentially decayed by a factor of0\.930\.93every2,0002\{,\}000iterations down to a floor of1×10−51\\times 10^\{\-5\}, and trained for60,00060\{,\}000iterations\.

During the warm\-up phase of FALM\-PINN \(Nwarm=5,000N\_\{\\mathrm\{warm\}\}=5\{,\}000\), all parameters are jointly trained by Adam at a fixed learning rate of2×10−32\\times 10^\{\-3\}\. The Fourier feature mapping usesD=800D=800features and bandwidthσ=1\\sigma=1\. The subsequent alternating optimization phase \(Nalter=55,000N\_\{\\mathrm\{alter\}\}=55\{,\}000\) updates the upper\-level parametersω\\omega, starting from a learning rate of1×10−41\\times 10^\{\-4\}and following the same exponential decay scheduler as the baselines\. Each lower\-level solve performsJ=2J=2LM iterations with the damping initialized atγ=1×10−5\\gamma=1\\times 10^\{\-5\}and decayed by a factor of0\.50\.5every2,0002\{,\}000iterations, until below1×10−61\\times 10^\{\-6\}\.

All baseline configurations are as follows\. RBA\[[31](https://arxiv.org/html/2608.05892#bib.bib31)\]usesηRBA=0\.001\\eta\_\{\\mathrm\{RBA\}\}=0\.001andγRBA=0\.999\\gamma\_\{\\mathrm\{RBA\}\}=0\.999\. CompleX\-PINN\[[40](https://arxiv.org/html/2608.05892#bib.bib40)\]adopts a single Cauchy layer of width200200\. PIKAN\[[43](https://arxiv.org/html/2608.05892#bib.bib43)\]uses a Kolmogorov–Arnold network of structure\[2,8,8,8,1\]\[2,8,8,8,1\]with grid size1010and spline order33\. SIREN\[[39](https://arxiv.org/html/2608.05892#bib.bib39)\]follows its default setting withω0=30\\omega\_\{0\}=30\.

Table[2](https://arxiv.org/html/2608.05892#S4.T2)and Fig\.[6](https://arxiv.org/html/2608.05892#S4.F6)summarize the numerical results for all methods\. During the joint warm\-up phase, FALM\-PINN exhibits a convergence trend comparable to that of the baseline PINNs\. After switching to the alternating phase, its relativeL2L^\{2\}error drops sharply and subsequently converges to4\.42×10−44\.42\\times 10^\{\-4\}\. This behavior indicates that the learned Fourier\-enhanced basis provides an effective representation, while the lower\-level LM updates rapidly refine the projection coefficients\. The results support the effectiveness of separating basis learning from coefficient fitting for nonlinear dispersive dynamics\.

Table 2:RelativeL2L^\{2\}errors on the 1D Korteweg–de Vries equation \(mean±\\pmstandard deviation over five independent runs\) for all methods\. The lowest error is shown in bold\.MethodRelativeL2L^\{2\}errorVanilla PINN4\.61×10−2±4\.03×10−34\.61\\times 10^\{\-2\}\\pm 4\.03\\times 10^\{\-3\}RBA7\.38×10−3±1\.32×10−37\.38\\times 10^\{\-3\}\\pm 1\.32\\times 10^\{\-3\}compleX\-PINN9\.63×10−3±1\.54×10−39\.63\\times 10^\{\-3\}\\pm 1\.54\\times 10^\{\-3\}PIKAN2\.52×10−2±7\.31×10−32\.52\\times 10^\{\-2\}\\pm 7\.31\\times 10^\{\-3\}SIREN2\.92×10−2±6\.18×10−32\.92\\times 10^\{\-2\}\\pm 6\.18\\times 10^\{\-3\}FALM\-PINN4\.42×𝟏𝟎−𝟒±1\.81×𝟏𝟎−𝟒\\mathbf\{4\.42\\times 10^\{\-4\}\\pm 1\.81\\times 10^\{\-4\}\}![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/kdv_comparison.png)Figure 6:RelativeL2L^\{2\}error trajectories on the 1D Korteweg–de Vries equation\. Curves denote the mean across five independent trials, and shaded bands indicate the min–max range\. The vertical dashed line marks the transition from the warm\-up phase \(Nwarm=5,000N\_\{\\mathrm\{warm\}\}=5\{,\}000\) to the alternating optimization phase\.The predicted solutions and corresponding point\-wise absolute errors are shown in Fig\.[7](https://arxiv.org/html/2608.05892#S4.F7)\. FALM\-PINN maintains uniformly lower errors across the spatio\-temporal domain and more accurately resolves the narrow oscillatory wave patterns\. These results demonstrate that the proposed method is effective in capturing the nonlinear dispersion structures\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/kdv_comparison_pred.png)Figure 7:Predictions \(left column\) and absolute error maps in log scale \(right column\) on the 1D Korteweg–de Vries equation for vanilla PINN, RBA, compleX\-PINN, PIKAN, SIREN, and FALM\-PINN respectively\.
### 4\.31D Heat equation with high\-Frequency solution

High\-frequency solutions are common failure modes for PINNs due to spectral bias\. To evaluate the ability of FALM\-PINN to capture oscillatory modes, we adopt the high\-frequency heat equation benchmark from\[[43](https://arxiv.org/html/2608.05892#bib.bib43)\]\. The equation is defined on the spatio\-temporal domainΩ=T×Ωs\\Omega=T\\times\\Omega\_\{s\}withT=\[0,1\]T=\[0,1\]andΩs=\[0,1\]\\Omega\_\{s\}=\[0,1\]:

ut=1\(F​π\)2​ux​x,\(t,x\)∈Ω,\\displaystyle u\_\{t\}=\\frac\{1\}\{\(F\\pi\)^\{2\}\}u\_\{xx\},\\quad\(t,x\)\\in\\Omega,\(45\)u​\(t,0\)=u​\(t,1\)=0,t∈T,\\displaystyle u\(t,0\)=u\(t,1\)=0,\\quad t\\in T,u​\(0,x\)=sin⁡\(F​π​x\),x∈Ωs,\\displaystyle u\(0,x\)=\\sin\(F\\pi x\),\\quad x\\in\\Omega\_\{s\},whereFFdenotes the spatial frequency\. The analytical solution isu​\(t,x\)=e−t​sin⁡\(F​π​x\)u\(t,x\)=e^\{\-t\}\\sin\(F\\pi x\), which combines a slowly decaying temporal component with a highly oscillatory spatial mode\. We setF=100F=100, which requires the network to recover high\-frequency structures and is known to be challenging for standard PINNs\. The exact solution is shown in Fig\.[8](https://arxiv.org/html/2608.05892#S4.F8)\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/heat1d.png)Figure 8:Exact solution of the 1D Heat equation\.All MLP\-based methods use a fully connected network with55hidden layers of width8080\. Training proceeds for50,00050\{,\}000iterations under Adam, with an initial learning rate of1\.5×10−31\.5\\times 10^\{\-3\}exponentially decayed by a factor of0\.80\.8every2,0002\{,\}000iterations down to1×10−51\\times 10^\{\-5\}\. The training set consists ofNf=15,000N\_\{f\}=15\{,\}000collocation points,Nb=200N\_\{b\}=200boundary points, andNi=200N\_\{i\}=200initial points\.

Since the heat equation is linear, the lower\-level problem of FALM\-PINN is convex inβ\\betaand admits a closed\-form one\-step update, so no warm\-up phase is required \(Nwarm=0N\_\{\\mathrm\{warm\}\}=0\)\. The Fourier feature mapping usesD=800D=800features with bandwidthσ=10\\sigma=10, and the damping is set toγ=1×10−7\\gamma=1\\times 10^\{\-7\}\.

For the baselines, we follow the high\-frequency tuning recommendations reported in their original papers\. RBA retainsηRBA=0\.001\\eta\_\{\\mathrm\{RBA\}\}=0\.001andγRBA=0\.999\\gamma\_\{\\mathrm\{RBA\}\}=0\.999\. CompleX\-PINN widens its single Cauchy layer to500500, while PIKAN employs a\[2,5,5,1\]\[2,5,5,1\]Kolmogorov–Arnold network with grid size150150and spline order33\. SIREN\[[39](https://arxiv.org/html/2608.05892#bib.bib39)\]employs sinusoidal activations withω0=100\\omega\_\{0\}=100, increased from the default value of3030to better capture the high\-frequency components present in the target solution\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/heat1d_comparison.png)Figure 9:Convergence trajectories on the 1D heat equation\. RelativeL2L^\{2\}error against training iterations; curves denote the mean across five independent trials, and shaded bands indicate the min–max range\.The training histories for all methods are shown in Fig\.[9](https://arxiv.org/html/2608.05892#S4.F9)\. None of the baselines resolves the high\-frequency mode within the50,00050\{,\}000iterations\. Vanilla PINN, RBA and compleX\-PINN fail in this approximation\. PIKAN and SIREN attain lower errors but still plateau at a high floor of about10−110^\{\-1\}\. In contrast, FALM\-PINN converges stably to an error of6\.8×10−46\.8\\times 10^\{\-4\}, more than two orders of magnitude below the best baseline\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/heat1d_pred_comparison.png)Figure 10:Predicted solutions \(left column\) and point\-wise absolute errors in log scale \(right column\) on the 1D heat equation for vanilla PINN, RBA, compleX\-PINN, PIKAN, SIREN, and FALM\-PINN, respectively\.As demonstrated in Fig\.[10](https://arxiv.org/html/2608.05892#S4.F10), vanilla PINN, RBA, and compleX\-PINN collapse to trivial, near\-zero predictions that miss the oscillation entirely\. PIKAN and SIREN show intermediate performance, capturing the general oscillatory pattern but with substantial point\-wise errors across the domain\. FALM\-PINN is the only method accurately recovering the high\-frequency mode with uniformly small error\. This consistent superiority confirms that FALM\-PINN can effectively mitigate spectral bias\.

To further assess the role of the Fourier\-enhanced basis in overcoming spectral bias, we conduct an ablation on the bandwidth hyperparameterσ\\sigma, which governs the frequency distribution of the induced kernel following Lemma[1](https://arxiv.org/html/2608.05892#Thmlemma1)\. We vary the spatial frequencyF∈\{10,100,150,200\}F\\in\\\{10,100,150,200\\\}against the bandwidthσ∈\{1,10,25,50\}\\sigma\\in\\\{1,10,25,50\\\}\. To ensure that the high\-frequency modes atF=200F=200are sufficiently sampled to satisfy the Nyquist\-Shannon criterion\[[60](https://arxiv.org/html/2608.05892#bib.bib60),[52](https://arxiv.org/html/2608.05892#bib.bib52)\], we double the number of boundary and initial points toNb=400N\_\{b\}=400andNi=400N\_\{i\}=400\. All remaining parameters are kept identical to the main experiment\.

The results of the ablation study are summarized in Table[3](https://arxiv.org/html/2608.05892#S4.T3); and Fig\.[11](https://arxiv.org/html/2608.05892#S4.F11)further visualizes the prediction slices at the terminal timet=1t=1across different spatial frequencies\. Both demonstrate that the bandwidth required for accurate approximation increases with the spatial frequencyFF\. For the low\-frequency case withF=10F=10, a narrow bandwidthσ=1\\sigma=1already suffices to span the spectral domain of the solution, and all tested bandwidths accurately reproduce the exact solution\. However, asFFgrows, only a sufficiently largeσ\\sigmaintroduces enough high\-frequency components into the Fourier\-enhanced features to capture the rapid oscillatory structure\. In the extreme regimeF=200F=200,σ=1\\sigma=1andσ=10\\sigma=10fail to recover meaningful predictions;σ=25\\sigma=25partially recovers the oscillation pattern but still exhibits noticeable approximation error; onlyσ=50\\sigma=50preserves an accurate approximation, achieving a relativeL2L^\{2\}error of3\.3×10−33\.3\\times 10^\{\-3\}\.

This trend is consistent with Lemma[1](https://arxiv.org/html/2608.05892#Thmlemma1), which indicates thatσ\\sigmashould scale with the dominant frequency of the target solution so that the Fourier\-enhanced features span the relevant spectral band\. Notably, even in the extremely high\-frequency regime withF=200F=200, a properly chosenσ\\sigmaenables FALM\-PINN to attain an accurate approximation, underscoring the role of the Fourier\-enhanced basis in mitigating spectral bias\.

Table 3:Ablation of the Fourier\-feature bandwidthσ\\sigmaacross spatial frequenciesFFon the 1D heat equation\. Each entry reports the relativeL2L^\{2\}error \(mean±\\pmstandard deviation\); the lowest error in each row is shown in bold\. A dash ’−\-’ denotes a configuration that failed to produce a meaningful prediction\.Bandwidthσ\\sigmaFrequencyFF1110102525505010108\.0×10−5±9\.9×10−68\.0\\times 10^\{\-5\}\\pm 9\.9\\times 10^\{\-6\}3\.6×𝟏𝟎−𝟓±7\.5×𝟏𝟎−𝟕\\mathbf\{3\.6\\times 10^\{\-5\}\\pm 7\.5\\times 10^\{\-7\}\}4\.7×10−5±2\.3×10−54\.7\\times 10^\{\-5\}\\pm 2\.3\\times 10^\{\-5\}5\.0×10−5±5\.0×10−65\.0\\times 10^\{\-5\}\\pm 5\.0\\times 10^\{\-6\}100100−\-6\.3×𝟏𝟎−𝟒±9\.5×𝟏𝟎−𝟓\\mathbf\{6\.3\\times 10^\{\-4\}\\pm 9\.5\\times 10^\{\-5\}\}1\.3×10−3±3\.5×10−41\.3\\times 10^\{\-3\}\\pm 3\.5\\times 10^\{\-4\}1\.0×10−3±1\.3×10−51\.0\\times 10^\{\-3\}\\pm 1\.3\\times 10^\{\-5\}150150−\-1\.9×10−1±1\.3×10−11\.9\\times 10^\{\-1\}\\pm 1\.3\\times 10^\{\-1\}2\.5×10−3±1\.1×10−32\.5\\times 10^\{\-3\}\\pm 1\.1\\times 10^\{\-3\}1\.8×𝟏𝟎−𝟑±3\.4×𝟏𝟎−𝟒\\mathbf\{1\.8\\times 10^\{\-3\}\\pm 3\.4\\times 10^\{\-4\}\}200200−\-−\-5\.9×10−1±8\.1×10−35\.9\\times 10^\{\-1\}\\pm 8\.1\\times 10^\{\-3\}3\.3×𝟏𝟎−𝟑±2\.4×𝟏𝟎−𝟒\\mathbf\{3\.3\\times 10^\{\-3\}\\pm 2\.4\\times 10^\{\-4\}\}

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/heat1d_freq_comparison_t1.png)Figure 11:Prediction slicesu​\(x,t=1\)u\(x,t=1\)of FALM\-PINN on the 1D heat equation for Fourier\-feature bandwidthsσ∈\{1,10,25,50\}\\sigma\\in\\\{1,10,25,50\\\}, at spatial frequenciesF∈\{200,150,100,10\}F\\in\\\{200,150,100,10\\\}\(top to bottom\)\. The black curve is the exact solution\.
### 4\.4Lid\-driven cavity flow

The lid\-driven cavity flow is a popular benchmark for steady incompressible Navier–Stokes equations with nonlinear convection and coupled velocity–pressure fields\. Following the setup in\[[61](https://arxiv.org/html/2608.05892#bib.bib61)\], we consider the system on the spatial domainΩs=\[0,1\]2\\Omega\_\{s\}=\[0,1\]^\{2\}:

\(𝐮⋅∇\)​𝐮\+∇p−1Re​Δ​𝐮\\displaystyle\(\\mathbf\{u\}\\cdot\\nabla\)\\mathbf\{u\}\+\\nabla p\-\\frac\{1\}\{\\mathrm\{Re\}\}\\Delta\\mathbf\{u\}=0,\\displaystyle=0,x\\displaystyle x∈Ωs,\\displaystyle\\in\\Omega\_\{s\},\(46\)∇⋅𝐮\\displaystyle\\nabla\\cdot\\mathbf\{u\}=0,\\displaystyle=0,x\\displaystyle x∈Ωs,\\displaystyle\\in\\Omega\_\{s\},where𝐮=\(u1,u2\)\\mathbf\{u\}=\(u\_\{1\},u\_\{2\}\)denotes the velocity field,ppdenotes the pressure, andRe=100\\mathrm\{Re\}=100is the Reynolds number\. Denoting the top boundary byΓ1\\Gamma\_\{1\}and the union of the left, right, and bottom boundaries byΓ2\\Gamma\_\{2\}, the flow is subject to

𝐮​\(x\)\\displaystyle\\mathbf\{u\}\(x\)=\(4​x​\(1−x\),0\),\\displaystyle=\\bigl\(4x\(1\-x\),0\\bigr\),x\\displaystyle x∈Γ1,\\displaystyle\\in\\Gamma\_\{1\},\(47\)𝐮​\(x\)\\displaystyle\\mathbf\{u\}\(x\)=\(0,0\),\\displaystyle=\(0,0\),x\\displaystyle x∈Γ2,\\displaystyle\\in\\Gamma\_\{2\},p​\(0,0\)\\displaystyle p\(0,0\)=0\.\\displaystyle=0\.
All MLP\-based methods use a fully connected network with33hidden layers, each with100100neurons\. We useNf=8,196N\_\{f\}=8\{,\}196collocation points andNb=256N\_\{b\}=256boundary points\. Training is performed using Adam for30,00030\{,\}000iterations, starting from an initial learning rate of1×10−31\\times 10^\{\-3\}, decayed by a factor of0\.990\.99every100100iterations\.

The warm\-up phase of FALM\-PINN lasts forNwarm=5,000N\_\{\\mathrm\{warm\}\}=5\{,\}000iterations via Adam at a fixed learning rate of1×10−31\\times 10^\{\-3\}\. The Fourier feature mapping usesD=800D=800features with bandwidthσ=1\\sigma=1\. During the alternating optimization phase \(Nalter=25,000N\_\{\\mathrm\{alter\}\}=25\{,\}000\), the upper\-level parametersω\\omegacontinue to be optimized using Adam with an initial learning rate of1×10−41\\times 10^\{\-4\}decayed by a factor of0\.990\.99every100100iterations down to a floor of1×10−51\\times 10^\{\-5\}\. Each lower\-level solve performsJ=2J=2LM iterations\. To accommodate the strong nonlinearity brought by the convection term, the LM damping is initialized atγ0=1×10−5\\gamma\_\{0\}=1\\times 10^\{\-5\}and decayed by a factor of0\.50\.5every2,0002\{,\}000iterations, until reaching1×10−61\\times 10^\{\-6\}\.

The remaining baseline configurations are as follows\. RBA applies residual\-based attention reweighting withηRBA=0\.001\\eta\_\{\\mathrm\{RBA\}\}=0\.001andγRBA=0\.999\\gamma\_\{\\mathrm\{RBA\}\}=0\.999\. CompleX\-PINN adopts a single Cauchy layer of width500500\. PIKAN uses a Kolmogorov–Arnold network of structure\[2,5,5,5,3\]\[2,5,5,5,3\]with grid size2020and spline order33\. SIREN employs sinusoidal activations withω0=30\\omega\_\{0\}=30as prescribed in the original paper\.

The relativeL2L^\{2\}error histories for all methods are presented in Fig\.[12](https://arxiv.org/html/2608.05892#S4.F12)\. During the warm\-up phase, the error reduction rate of FALM\-PINN remains on the same scale as the baselines\. Upon transition to the alternating optimization phase, the error decreases by nearly one order of magnitude and continues to decay smoothly, which is consistent with the two\-stage behavior observed in the 2D Klein–Gordon experiment\. Table[4](https://arxiv.org/html/2608.05892#S4.T4)further reports the component\-wise relativeL2L^\{2\}errors foru1u\_\{1\},u2u\_\{2\}, andpp\. FALM\-PINN achieves the lowest error in all three components and reduces the total error by about45%45\\%relative to SIREN, the most accurate baseline in this experiment\.

Table 4:RelativeL2L^\{2\}error \(mean±\\pmstandard deviation\) for each solution component \(u1u\_\{1\},u2u\_\{2\},pp\) on the 2D lid\-driven cavity flow\. The lowest error in each column is shown in bold\.Methodu1u\_\{1\}u2u\_\{2\}ppTotalVanilla PINN2\.54×10−2±5\.63×10−32\.54\\times 10^\{\-2\}\\pm 5\.63\\times 10^\{\-3\}4\.81×10−2±1\.21×10−24\.81\\times 10^\{\-2\}\\pm 1\.21\\times 10^\{\-2\}8\.02×10−2±2\.43×10−28\.02\\times 10^\{\-2\}\\pm 2\.43\\times 10^\{\-2\}4\.48×10−2±7\.75×10−34\.48\\times 10^\{\-2\}\\pm 7\.75\\times 10^\{\-3\}RBA1\.72×10−2±1\.53×10−31\.72\\times 10^\{\-2\}\\pm 1\.53\\times 10^\{\-3\}2\.55×10−2±1\.83×10−32\.55\\times 10^\{\-2\}\\pm 1\.83\\times 10^\{\-3\}3\.69×10−2±4\.93×10−33\.69\\times 10^\{\-2\}\\pm 4\.93\\times 10^\{\-3\}2\.23×10−2±1\.89×10−32\.23\\times 10^\{\-2\}\\pm 1\.89\\times 10^\{\-3\}compleX\-PINN1\.28×10−2±6\.94×10−41\.28\\times 10^\{\-2\}\\pm 6\.94\\times 10^\{\-4\}2\.22×10−2±1\.71×10−32\.22\\times 10^\{\-2\}\\pm 1\.71\\times 10^\{\-3\}3\.71×10−2±2\.44×10−33\.71\\times 10^\{\-2\}\\pm 2\.44\\times 10^\{\-3\}1\.96×10−2±1\.16×10−31\.96\\times 10^\{\-2\}\\pm 1\.16\\times 10^\{\-3\}PIKAN4\.99×10−3±8\.37×10−44\.99\\times 10^\{\-3\}\\pm 8\.37\\times 10^\{\-4\}7\.72×10−3±1\.08×10−37\.72\\times 10^\{\-3\}\\pm 1\.08\\times 10^\{\-3\}1\.57×10−2±2\.53×10−31\.57\\times 10^\{\-2\}\\pm 2\.53\\times 10^\{\-3\}7\.82×10−3±1\.43×10−37\.82\\times 10^\{\-3\}\\pm 1\.43\\times 10^\{\-3\}SIREN3\.11×10−3±5\.49×10−43\.11\\times 10^\{\-3\}\\pm 5\.49\\times 10^\{\-4\}4\.28×10−3±1\.33×10−34\.28\\times 10^\{\-3\}\\pm 1\.33\\times 10^\{\-3\}1\.31×10−2±1\.30×10−31\.31\\times 10^\{\-2\}\\pm 1\.30\\times 10^\{\-3\}5\.07×10−3±4\.59×10−45\.07\\times 10^\{\-3\}\\pm 4\.59\\times 10^\{\-4\}FALM\-PINN1\.11×𝟏𝟎−𝟑±3\.38×𝟏𝟎−𝟒\\mathbf\{1\.11\\times 10^\{\-3\}\\pm 3\.38\\times 10^\{\-4\}\}2\.06×𝟏𝟎−𝟑±1\.87×𝟏𝟎−𝟒\\mathbf\{2\.06\\times 10^\{\-3\}\\pm 1\.87\\times 10^\{\-4\}\}8\.04×𝟏𝟎−𝟑±6\.41×𝟏𝟎−𝟒\\mathbf\{8\.04\\times 10^\{\-3\}\\pm 6\.41\\times 10^\{\-4\}\}2\.78×𝟏𝟎−𝟑±1\.29×𝟏𝟎−𝟒\\mathbf\{2\.78\\times 10^\{\-3\}\\pm 1\.29\\times 10^\{\-4\}\}

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/ns2d_comparison.png)Figure 12:RelativeL2L^\{2\}error trajectories on the 2D lid\-driven cavity flow\. Curves denote the mean across five independent trials; shaded bands indicate the min–max range\. The vertical dashed line marks the transition from the warm\-up phase \(Nwarm=5,000N\_\{\\mathrm\{warm\}\}=5\{,\}000\) to the alternating optimization phase\.We further examine the point\-wise absolute errors of the predicted velocity componentsu1u\_\{1\},u2u\_\{2\}, and pressureppin Fig\.[14](https://arxiv.org/html/2608.05892#S4.F14), computed against the reference solution shown in Fig\.[13](https://arxiv.org/html/2608.05892#S4.F13)\. FALM\-PINN achieves smaller errors across all three fields, indicating that its improvement is reflected not only in the relativeL2L^\{2\}error but also in the spatial error distribution\. These results show that FALM\-PINN accurately reproduces the coupled flow structure and remains effective for approximating nonlinear multi\-component PDE systems\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/ns2d.png)Figure 13:Reference solution of the 2D lid\-driven cavity flow\.![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/ns2d_solution_comparison.png)Figure 14:Point\-wise absolute error maps \(log scale\) on the 2D lid\-driven cavity flow\. Columns: horizontal velocityu1u\_\{1\}, vertical velocityu2u\_\{2\}, and pressurepp\. Rows: vanilla PINN, RBA, compleX\-PINN, PIKAN, SIREN, and FALM\-PINN\.
### 4\.51D viscous Burgers equation

The viscous Burgers equation is a nonlinear convection–diffusion benchmark that develops shock\-like steep\-gradient structures under small viscosity\. Following the formulation in\[[12](https://arxiv.org/html/2608.05892#bib.bib12)\], we consider the equation on the spatio\-temporal domainΩ=T×Ωs\\Omega=T\\times\\Omega\_\{s\}withΩs=\[−1,1\]\\Omega\_\{s\}=\[\-1,1\]andT=\[0,1\]T=\[0,1\]:

ut\+u​ux−0\.01π​ux​x=0,\(t,x\)∈Ω,\\displaystyle u\_\{t\}\+uu\_\{x\}\-\\frac\{0\.01\}\{\\pi\}u\_\{xx\}=0,\\quad\(t,x\)\\in\\Omega,\(48\)u​\(0,x\)=−sin⁡\(π​x\),x∈Ωs,\\displaystyle u\(0,x\)=\-\\sin\(\\pi x\),\\quad x\\in\\Omega\_\{s\},u​\(t,−1\)=u​\(t,1\)=0,t∈T\.\\displaystyle u\(t,\-1\)=u\(t,1\)=0,\\quad t\\in T\.All MLP\-based methods share the same network architecture of88hidden layers of width2020\. We useNf=10,000N\_\{f\}=10\{,\}000collocation points andNi\+Nb=100N\_\{i\}\+N\_\{b\}=100combined initial and boundary points\. The baselines are optimized by Adam with an initial learning rate of5×10−35\\times 10^\{\-3\}, exponentially decayed by a factor of0\.70\.7every1,0001\{,\}000iterations down to a floor of1×10−51\\times 10^\{\-5\}, and trained for40,00040\{,\}000iterations for all methods to fully converge\.

During the warm\-up phase of FALM\-PINN \(Nwarm=10,000N\_\{\\mathrm\{warm\}\}=10\{,\}000\), all parameters are jointly trained by Adam at a fixed learning rate of1×10−31\\times 10^\{\-3\}\. The Fourier feature mapping usesD=800D=800features and bandwidthσ=1\\sigma=1\. The subsequent alternating optimization phase \(Nalter=30,000N\_\{\\mathrm\{alter\}\}=30\{,\}000\) then updates the upper\-level parametersω\\omega, starting from a learning rate of1×10−41\\times 10^\{\-4\}and following the same exponential decay scheduler as the baselines\. Each lower\-level solve performsJ=1J=1LM iteration with the damping initialized atγ=1×10−5\\gamma=1\\times 10^\{\-5\}and decayed by a factor of0\.50\.5for every10001000iterations, until below1×10−61\\times 10^\{\-6\}\.

For the baselines, RBA\[[31](https://arxiv.org/html/2608.05892#bib.bib31)\]uses residual\-based attention reweighting withηRBA=0\.001\\eta\_\{\\mathrm\{RBA\}\}=0\.001andγRBA=0\.999\\gamma\_\{\\mathrm\{RBA\}\}=0\.999\. CompleX\-PINN\[[40](https://arxiv.org/html/2608.05892#bib.bib40)\]adopts a single Cauchy layer of width200200\. PIKAN\[[43](https://arxiv.org/html/2608.05892#bib.bib43)\]uses a Kolmogorov–Arnold network of structure\[2,5,5,5,1\]\[2,5,5,5,1\]with grid size1515and spline order33\. SIREN\[[39](https://arxiv.org/html/2608.05892#bib.bib39)\]is used with its default configuration withω0=30\\omega\_\{0\}=30\.

Table[5](https://arxiv.org/html/2608.05892#S4.T5)and Fig\.[15](https://arxiv.org/html/2608.05892#S4.F15)summarize the final results and the convergence behavior of all methods\. During the warm\-up phase, FALM\-PINN follows a decay trend comparable to the baselines\. After switching to the alternating optimization phase, however, the error undergoes a sharp decline and reaches a final relativeL2L^\{2\}error of5\.10×10−55\.10\\times 10^\{\-5\}, while all baseline methods remain above1×10−31\\times 10^\{\-3\}\. This two\-stage pattern echoes the behavior observed in the previous experiments and further supports the effectiveness of the framework\. However, the training time and memory cost of FALM\-PINN are higher than those of the baselines on this problem\. This is because FALM\-PINN relies on the residual matrix, whose dimension grows with the number of collocation points\. Consequently, its computational and memory costs increase more rapidly than the baselines as the dataset size grows\.

Table 5:1D viscous Burgers equation: relativeL2L^\{2\}errors \(mean±\\pmstd over five independent runs\), GPU time per10001000iterations, and peak GPU memory for all methods\. The lowest value in each column is shown in bold\.MethodRelativeL2L^\{2\}errorTime \(s\)Memory \(GB\)Vanilla PINN2\.16×10−2±5\.99×10−32\.16\\times 10^\{\-2\}\\pm 5\.99\\times 10^\{\-3\}3\.29\\mathbf\{3\.29\}0\.08\\mathbf\{0\.08\}RBA3\.03×10−3±9\.16×10−43\.03\\times 10^\{\-3\}\\pm 9\.16\\times 10^\{\-4\}3\.953\.950\.090\.09compleX\-PINN3\.53×10−2±4\.47×10−33\.53\\times 10^\{\-2\}\\pm 4\.47\\times 10^\{\-3\}4\.034\.030\.440\.44PIKAN2\.13×10−2±7\.69×10−32\.13\\times 10^\{\-2\}\\pm 7\.69\\times 10^\{\-3\}15\.5915\.590\.430\.43SIREN1\.53×10−2±5\.65×10−31\.53\\times 10^\{\-2\}\\pm 5\.65\\times 10^\{\-3\}4\.124\.120\.100\.10FALM\-PINN5\.10×𝟏𝟎−𝟓±1\.22×𝟏𝟎−𝟓\\mathbf\{5\.10\\times 10^\{\-5\}\\pm 1\.22\\times 10^\{\-5\}\}14\.3514\.351\.251\.25![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/comparison_burgers1d_plot.png)Figure 15:RelativeL2L^\{2\}error trajectories on the 1D viscous Burgers equation\. Curves denote the mean across five independent trials, and shaded bands indicate the min–max range\.The reference solution and the corresponding point\-wise absolute errors are presented in Figs\.[16](https://arxiv.org/html/2608.05892#S4.F16)and[17](https://arxiv.org/html/2608.05892#S4.F17)\. FALM\-PINN maintains lower errors across the domain, including the steep\-gradient region that develops nearx=0x=0\. These results indicate that the improvement is not limited to the global relativeL2L^\{2\}error but also appears in the localized region where the nonlinear convection term makes the solution difficult to approximate\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/burgers1d.png)Figure 16:Reference solution of the 1D viscous Burgers equation\.![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/burgers1d_solution_comparison.png)Figure 17:Predictions \(left column\) and error maps in log scale \(right column\) on the 1D viscous Burgers equation for vanilla PINN, RBA, compleX\-PINN, PIKAN, SIREN and FALM\-PINN respectively\.To better understand the behavior of all predictions near the shock, we present temporal prediction slices for all methods at fixed spatial locationsx∈\{−0\.1,0,0\.1\}x\\in\\\{\-0\.1,0,0\.1\\\}against the reference solution in Fig\.[18](https://arxiv.org/html/2608.05892#S4.F18)\. FALM\-PINN closely tracks the reference at all three locations throughout the temporal domain\. In contrast, all baselines reproduce only the overall profile and exhibit visible deviations in the critical region neart=0\.5t=0\.5\. These slices further validate the robustness of FALM\-PINN in resolving nonlinear PDEs with shock structure\.

![Refer to caption](https://arxiv.org/html/2608.05892v1/Figure/burgers1d_x_slices.png)Figure 18:Prediction slices over time at fixed spatial positionsx∈\{−0\.1,0,0\.1\}x\\in\\\{\-0\.1,0,0\.1\\\}on the 1D viscous Burgers equation, comparing all methods against the exact solution\.
### 4\.6Comparison with IFeF\-PINN

Recall that IFeF\-PINN is designed primarily for scalar, linear PDEs, where the lower\-level problem is convex and admits a closed\-form solution\[[52](https://arxiv.org/html/2608.05892#bib.bib52)\]\. For nonlinear PDEs, this problem becomes nonconvex and no longer has a closed\-form solution; IFeF\-PINN therefore updates the coefficients to a local minimizer everyNlowerN\_\{\\text\{lower\}\}epochs using L\-BFGS\. As a generic optimizer, it captures the overall curvature of the loss but ignores the least\-squares structure of the residual\. To assess the effectiveness of the lower\-level LM solver of FALM\-PINN, we compare the two methods on the three nonlinear benchmarks \(2D Klein–Gordon, 1D Korteweg–de Vries, and 1D viscous Burgers\), keeping the network architecture, the Fourier\-feature mapping, the training set, and all training hyperparameters identical\.

Table[6](https://arxiv.org/html/2608.05892#S4.T6)reports the final relativeL2L^\{2\}errors\. FALM\-PINN is one to two orders of magnitude more accurate than IFeF\-PINN on every nonlinear benchmark\. The L\-BFGS lower\-level of IFeF\-PINN does reduce the loss, but it converges slowly and stalls at a substantially higher error, because it does not exploit the least\-squares geometry that the LM update captures through the Gauss–Newton Hessian and adaptive damping\. These results show that, for the nonlinear decoupled problem, the structure\-exploiting LM solve is markedly more accurate than a generic gradient\-based optimizer, even though both share the same Fourier\-enhanced basis\.

Table 6:Comparison of IFeF\-PINN and FALM\-PINN on three scalar nonlinear PDE benchmarks\. All settings other than the lower\-level solver are kept identical\. Entries report the final relativeL2L^\{2\}error \(mean±\\pmstd over five independent runs\)\.BenchmarkIFeF\-PINNFALM\-PINN2D Klein–Gordon1\.28×10−3±1\.87×10−41\.28\\times 10^\{\-3\}\\pm 1\.87\\times 10^\{\-4\}3\.50×𝟏𝟎−𝟓±4\.88×𝟏𝟎−𝟔\\mathbf\{3\.50\\times 10^\{\-5\}\\pm 4\.88\\times 10^\{\-6\}\}1D Korteweg–de Vries1\.32×10−2±4\.28×10−31\.32\\times 10^\{\-2\}\\pm 4\.28\\times 10^\{\-3\}4\.42×𝟏𝟎−𝟒±1\.81×𝟏𝟎−𝟒\\mathbf\{4\.42\\times 10^\{\-4\}\\pm 1\.81\\times 10^\{\-4\}\}1D viscous Burgers2\.46×10−3±6\.71×10−42\.46\\times 10^\{\-3\}\\pm 6\.71\\times 10^\{\-4\}5\.10×𝟏𝟎−𝟓±1\.22×𝟏𝟎−𝟓\\mathbf\{5\.10\\times 10^\{\-5\}\\pm 1\.22\\times 10^\{\-5\}\}

## 5Conclusion and future work

In this work, we proposed FALM\-PINN, an alternating training framework that decouples Fourier\-enhanced basis learning from projection\-coefficient fitting\. The upper\-level problem learns an adaptive latent basis, while the lower\-level problem exploits the nonlinear least\-squares structure of the physics residual through damped Levenberg–Marquardt updates\. We establish global convergence for this formulation, overcoming the limitation of conventional PINNs\. Empirically, FALM\-PINN significantly improves accuracy compared to state\-of\-the\-art methods across a variety of high\-frequency, nonlinear, and coupled benchmarks, achieving substantially lower relativeL2L^\{2\}errors than the considered baselines\.

A promising direction for future work is to develop adaptive strategies for selecting the Fourier bandwidth, thereby constructing a more compact basis without sacrificing approximation accuracy\. Combining FALM\-PINN with adaptive resampling methods may further reduce the computational cost and memory consumption of the method\. These developments would improve the scalability of the framework to more complex and high\-dimensional PDEs\. In parallel, integrating FALM\-PINN with advanced neural architectures or operator\-learning frameworks offers further opportunities to enhance its expressiveness and generalization capability\. Such extensions may lead to more accurate and reliable physics\-informed solvers and help narrow the performance gap between learning\-based approaches and classical numerical methods\.

## Appendix AWorked example: scalar conservation laws

To illustrate the assembly of the residual vector and its Jacobian in the lower\-level solver, we instantiate the framework of Section[3\.2](https://arxiv.org/html/2608.05892#S3.SS2)on the scalar conservation law:

ut\+\[P​\(u\)\]x=0,\(t,x\)∈\(0,T\]×\(0,1\),u\_\{t\}\+\\bigl\[P\(u\)\\bigr\]\_\{x\}=0,\\quad\(t,x\)\\in\(0,T\]\\times\(0,1\),\(49\)subject to Dirichlet boundary conditionsu​\(t,0\)=u​\(t,1\)=0u\(t,0\)=u\(t,1\)=0and initial conditionu​\(0,x\)=u0​\(x\)u\(0,x\)=u\_\{0\}\(x\), whereP:ℝ→ℝP:\\mathbb\{R\}\\to\\mathbb\{R\}denotes a smooth nonlinear flux function\. This scalar form covers a wide range of one\-dimensional nonlinear conservation laws\. Representative examples includeP​\(u\)=12​u2P\(u\)=\\tfrac\{1\}\{2\}u^\{2\}, which leads to the inviscid Burgers’ equation, andP​\(u\)=u​e−uP\(u\)=ue^\{\-u\}, commonly used in traffic flow models\.

For the scalar conservation law \([49](https://arxiv.org/html/2608.05892#A1.E49)\) withm=M=1m=M=1, the coefficient vector reduces toβ∈ℝ2​D\\beta\\in\\mathbb\{R\}^\{2D\}and the approximate solution at thejj\-th step isuj​\(x\)=ψD​\(x\)⊤​β\(j\)u\_\{j\}\(x\)=\\psi\_\{D\}\(x\)^\{\\top\}\\beta^\{\(j\)\}\. By the chain rule\[P​\(u\)\]x=P′​\(u\)​ux\[P\(u\)\]\_\{x\}=P^\{\\prime\}\(u\)u\_\{x\}, the nonlinear PDE operator at a collocation pointxfix\_\{f\}^\{i\}gives:

𝒩​\[uω,β\]​\(xfi\)\\displaystyle\\mathscr\{N\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{i\}\)=\(∂tψDi\)⊤​β\(j\)\+P′​\(uji\)​\(∂xψDi\)⊤​β\(j\),\\displaystyle=\(\\partial\_\{t\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\\beta^\{\(j\)\}\+P^\{\\prime\}\(u\_\{j\}^\{i\}\)\(\\partial\_\{x\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\\beta^\{\(j\)\},\(50\)𝑴f​\(β\)\\displaystyle\\bm\{M\}\_\{f\}\(\\beta\)=\[𝒩​\[uω,β\]​\(xf1\),…,𝒩​\[uω,β\]​\(xfNf\)\]⊤,\\displaystyle=\\big\[\\mathscr\{N\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{1\}\),\\dots,\\mathscr\{N\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{N\_\{f\}\}\)\\big\]^\{\\top\},where𝑴f​\(β\)∈ℝNf\\bm\{M\}\_\{f\}\(\\beta\)\\in\\mathbb\{R\}^\{N\_\{f\}\}and we abbreviateψDi=ψD​\(xfi\)\\psi\_\{D\}^\{i\}=\\psi\_\{D\}\(x\_\{f\}^\{i\}\),uji=ψDi⊤​β\(j\)u\_\{j\}^\{i\}=\{\\psi\_\{D\}^\{i\}\}^\{\\top\}\\beta^\{\(j\)\}, and∂tψDi\\partial\_\{t\}\\psi\_\{D\}^\{i\},∂xψDi\\partial\_\{x\}\\psi\_\{D\}^\{i\}denotes the partial derivatives of the feature map\. Differentiating𝒩​\[uω,β\]​\(xfi\)\\mathscr\{N\}\[u\_\{\\omega,\\beta\}\]\(x\_\{f\}^\{i\}\)with respect toβ\\betayields the Jacobian matrix∇β𝑴f​\(β\(j\)\)∈ℝNf×2​D\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)\\in\\mathbb\{R\}^\{N\_\{f\}\\times 2D\}:

∇β𝑴f​\(β\(j\)\)=\[\(∂tψD1\)⊤\+P′​\(uj1\)​\(∂xψD1\)⊤\+P′′​\(uj1\)​\(∂xψD1\)⊤​β\(j\)​\(ψD1\)⊤⋮\(∂tψDNf\)⊤\+P′​\(ujNf\)​\(∂xψDNf\)⊤\+P′′​\(ujNf\)​\(∂xψDNf\)⊤​β\(j\)​\(ψDNf\)⊤\]\.\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)=\\begin\{bmatrix\}\(\\partial\_\{t\}\\psi\_\{D\}^\{1\}\)^\{\\top\}\+P^\{\\prime\}\(u\_\{j\}^\{1\}\)\(\\partial\_\{x\}\\psi\_\{D\}^\{1\}\)^\{\\top\}\+P^\{\\prime\\prime\}\(u\_\{j\}^\{1\}\)\(\\partial\_\{x\}\\psi\_\{D\}^\{1\}\)^\{\\top\}\\beta^\{\(j\)\}\(\\psi\_\{D\}^\{1\}\)^\{\\top\}\\\\ \\vdots\\\\ \(\\partial\_\{t\}\\psi\_\{D\}^\{N\_\{f\}\}\)^\{\\top\}\+P^\{\\prime\}\(u\_\{j\}^\{N\_\{f\}\}\)\(\\partial\_\{x\}\\psi\_\{D\}^\{N\_\{f\}\}\)^\{\\top\}\+P^\{\\prime\\prime\}\(u\_\{j\}^\{N\_\{f\}\}\)\(\\partial\_\{x\}\\psi\_\{D\}^\{N\_\{f\}\}\)^\{\\top\}\\beta^\{\(j\)\}\(\\psi\_\{D\}^\{N\_\{f\}\}\)^\{\\top\}\\end\{bmatrix\}\.\(51\)For the inviscid Burgers’ fluxP​\(u\)=12​u2P\(u\)=\\tfrac\{1\}\{2\}u^\{2\}, we haveP′​\(uji\)=ujiP^\{\\prime\}\(u\_\{j\}^\{i\}\)=u\_\{j\}^\{i\}andP′′​\(uji\)=1P^\{\\prime\\prime\}\(u\_\{j\}^\{i\}\)=1\. Consequently, each row of the physics Jacobian∇β𝑴f​\(β\(j\)\)\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)simplifies to:

\(∂tψDi\)⊤\+uji​\(∂xψDi\)⊤\+\(\(∂xψDi\)⊤​β\(j\)\)​\(ψDi\)⊤∈ℝ1×2​D\.\(\\partial\_\{t\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\+u\_\{j\}^\{i\}\(\\partial\_\{x\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\+\\bigl\(\(\\partial\_\{x\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\\beta^\{\(j\)\}\\bigr\)\(\\psi\_\{D\}^\{i\}\)^\{\\top\}\\in\\mathbb\{R\}^\{1\\times 2D\}\.\(52\)For the traffic flow modelP​\(u\)=u​e−uP\(u\)=ue^\{\-u\}, we haveP′​\(uji\)=\(1−uji\)​e−ujiP^\{\\prime\}\(u\_\{j\}^\{i\}\)=\(1\-u\_\{j\}^\{i\}\)e^\{\-u\_\{j\}^\{i\}\}andP′′​\(uji\)=\(uji−2\)​e−ujiP^\{\\prime\\prime\}\(u\_\{j\}^\{i\}\)=\(u\_\{j\}^\{i\}\-2\)e^\{\-u\_\{j\}^\{i\}\}\. Accordingly, each row of∇β𝑴f​\(β\(j\)\)\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)becomes:

\(∂tψDi\)⊤\+\(1−uji\)​e−uji​\(∂xψDi\)⊤\+\(uji−2\)​e−uji​\(\(∂xψDi\)⊤​β\(j\)\)​\(ψDi\)⊤∈ℝ1×2​D\.\(\\partial\_\{t\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\+\(1\-u\_\{j\}^\{i\}\)e^\{\-u\_\{j\}^\{i\}\}\(\\partial\_\{x\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\+\(u\_\{j\}^\{i\}\-2\)e^\{\-u\_\{j\}^\{i\}\}\\bigl\(\(\\partial\_\{x\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\\beta^\{\(j\)\}\\bigr\)\(\\psi\_\{D\}^\{i\}\)^\{\\top\}\\in\\mathbb\{R\}^\{1\\times 2D\}\.\(53\)For a linear fluxP​\(u\)=a​uP\(u\)=au,P′​\(uji\)=aP^\{\\prime\}\(u\_\{j\}^\{i\}\)=aandP′′​\(uji\)=0P^\{\\prime\\prime\}\(u\_\{j\}^\{i\}\)=0, so each row becomes\(∂tψDi\)⊤\+a​\(∂xψDi\)⊤\(\\partial\_\{t\}\\psi\_\{D\}^\{i\}\)^\{\\top\}\+a\(\\partial\_\{x\}\\psi\_\{D\}^\{i\}\)^\{\\top\}, which is independent ofβ\(j\)\\beta^\{\(j\)\}\.

With the Jacobian assembled, the updateβ\(j\+1\)=β\(j\)\+Δ​β\(j\)\\beta^\{\(j\+1\)\}=\\beta^\{\(j\)\}\+\\Delta\\beta^\{\(j\)\}is obtained by solving the normal equation \([18](https://arxiv.org/html/2608.05892#S3.E18)\), where the Gram matrix expands as:

𝐉j⊤​𝐉j=1Nb​𝑴b⊤​𝑴b\+λNf​\(∇β𝑴f​\(β\(j\)\)\)⊤​∇β𝑴f​\(β\(j\)\)∈ℝ2​D×2​D\.\\mathbf\{J\}\_\{j\}^\{\\top\}\\mathbf\{J\}\_\{j\}=\\frac\{1\}\{N\_\{b\}\}\\bm\{M\}\_\{b\}^\{\\top\}\\bm\{M\}\_\{b\}\+\\frac\{\\lambda\}\{N\_\{f\}\}\\bigl\(\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)\\bigr\)^\{\\top\}\\nabla\_\{\\beta\}\\bm\{M\}\_\{f\}\(\\beta^\{\(j\)\}\)\\in\\mathbb\{R\}^\{2D\\times 2D\}\.\(54\)The boundary Gram matrix𝑴b⊤​𝑴b\\bm\{M\}\_\{b\}^\{\\top\}\\bm\{M\}\_\{b\}is constant, while the physics Gram matrix is recomputed at each iteration\.

## References

- \[1\]R\. J\. LeVeque, Finite difference methods for ordinary and partial differential equations: steady\-state and time\-dependent problems, SIAM, 2007\.
- \[2\]O\. C\. Zienkiewicz, R\. L\. Taylor, P\. Nithiarasu, J\. Zhu, The finite element method, Vol\. 3, Elsevier, 1977\.
- \[3\]P\. G\. Ciarlet, The finite element method for elliptic problems, SIAM, 2002\.
- \[4\]S\. C\. Brenner, L\. R\. Scott, The mathematical theory of finite element methods, Springer, 2008\.
- \[5\]J\. P\. Boyd, Chebyshev and Fourier spectral methods, Courier Corporation, 2001\.
- \[6\]C\. Canuto, M\. Y\. Hussaini, A\. Quarteroni, T\. A\. Zang, Spectral methods, Vol\. 285, Springer, 2006\.
- \[7\]J\. Shen, T\. Tang, L\.\-L\. Wang, Spectral methods: algorithms, analysis and applications, Vol\. 41, Springer Science & Business Media, 2011\.
- \[8\]E\. Weinan, B\. Engquist, The heterognous multiscale methods, Communications in Mathematical Sciences 1 \(1\) \(2003\) 87–132\.
- \[9\]A\. Abdulle, E\. Weinan, B\. Engquist, E\. Vanden\-Eijnden, The heterogeneous multiscale method, Acta Numerica 21 \(2012\) 1–87\.
- \[10\]Y\. Efendiev, T\. Y\. Hou, Multiscale finite element methods: theory and applications, Springer Science & Business Media, 2009\.
- \[11\]A\. Iserles, S\. P\. Nørsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461 \(2057\) \(2005\) 1383–1399\.
- \[12\]M\. Raissi, P\. Perdikaris, G\. E\. Karniadakis, Physics\-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 \(2019\) 686–707\.
- \[13\]G\. E\. Karniadakis, I\. G\. Kevrekidis, L\. Lu, P\. Perdikaris, S\. Wang, L\. Yang, Physics\-informed machine learning, Nature Reviews Physics \(2021\)\.
- \[14\]Z\. Hu, K\. Shukla, G\. E\. Karniadakis, K\. Kawaguchi, Tackling the curse of dimensionality with physics\-informed neural networks, Neural Networks 176 \(2024\) 106369\.
- \[15\]F\. S\. Costabal, S\. Pezzuto, P\. Perdikaris,δ\\delta\-PINNs: Physics\-informed neural networks on complex geometries, Engineering Applications of Artificial Intelligence 127 \(2024\) 107324\.
- \[16\]M\. Raissi, A\. Yazdani, G\. E\. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367 \(6481\) \(2020\) 1026–1030\.
- \[17\]X\. Jin, S\. Cai, H\. Li, G\. E\. Karniadakis, NSFnets \(Navier\-Stokes flow nets\): Physics\-informed neural networks for the incompressible Navier\-Stokes equations, Journal of Computational Physics 426 \(2021\) 109951\.
- \[18\]Y\. Chen, L\. Lu, G\. E\. Karniadakis, L\. Dal Negro, Physics\-informed neural networks for inverse problems in nano\-optics and metamaterials, Optics express 28 \(8\) \(2020\) 11618–11633\.
- \[19\]L\. Yang, X\. Meng, G\. E\. Karniadakis, B\-pinns: Bayesian physics\-informed neural networks for forward and inverse pde problems with noisy data, Journal of Computational Physics 425 \(2021\) 109913\.
- \[20\]J\.\-H\. Bastek, W\. Sun, D\. Kochmann, Physics\-informed diffusion models, in: International Conference on Learning Representations, Vol\. 2025, 2025, pp\. 3360–3385\.
- \[21\]Z\. Wang, A\. Harting, M\. Barreau, M\. M\. Zavlanos, K\. H\. Johansson, Source\-guided flow matching, arXiv preprint arXiv:2508\.14807 \(2025\)\.
- \[22\]D\. Shu, Z\. Li, A\. B\. Farimani, A physics\-informed diffusion model for high\-fidelity flow field reconstruction, Journal of Computational Physics 478 \(2023\) 111972\.
- \[23\]Z\. Li, N\. B\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, A\. Anandkumar, Fourier neural operator for parametric partial differential equations, in: International Conference on Learning Representations, 2021\.
- \[24\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, G\. E\. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence 3 \(3\) \(2021\) 218–229\.
- \[25\]N\. Rahaman, A\. Baratin, D\. Arpit, F\. Draxler, M\. Lin, F\. Hamprecht, Y\. Bengio, A\. Courville, On the spectral bias of neural networks, in: International Conference on Machine Learning, PMLR, 2019, pp\. 5301–5310\.
- \[26\]Z\.\-Q\. J\. Xu, L\. Zhang, W\. Cai, On understanding and overcoming spectral biases of deep neural network learning methods for solving PDEs, Journal of Computational Physics \(2025\) 113905\.
- \[27\]A\. Krishnapriyan, A\. Gholami, S\. Zhe, R\. Kirby, M\. W\. Mahoney, Characterizing possible failure modes in physics\-informed neural networks, Advances in Neural Information Processing Systems 34 \(2021\) 26548–26560\.
- \[28\]S\. Wang, S\. Sankaran, P\. Perdikaris, Respecting causality for training physics\-informed neural networks, Computer Methods in Applied Mechanics and Engineering 421 \(2024\) 116813\.
- \[29\]S\. Wang, X\. Yu, P\. Perdikaris, When and why PINNs fail to train: A neural tangent kernel perspective, Journal of Computational Physics 449 \(2022\) 110768\.
- \[30\]Y\. Song, H\. Wang, H\. Yang, M\. L\. Taccari, X\. Chen, Loss\-attentional physics\-informed neural networks, Journal of Computational Physics 501 \(2024\) 112781\.
- \[31\]S\. J\. Anagnostopoulos, J\. D\. Toscano, N\. Stergiopulos, G\. E\. Karniadakis, Residual\-based attention in physics\-informed neural networks, Computer Methods in Applied Mechanics and Engineering 421 \(2024\) 116805\.
- \[32\]C\. Si, M\. Yan, Convolution\-weighting method for the physics\-informed neural network: A primal\-dual optimization perspective, Journal of Computational Physics 555 \(2026\) 113911\.
- \[33\]C\. Zhao, X\. Xie, W\. Chen, Casual attention: Adaptive enforcement of causality in physics\-informed neural networks, Journal of Computational Physics \(2026\) 115071\.
- \[34\]C\. Wu, M\. Zhu, Q\. Tan, Y\. Kartha, L\. Lu, A comprehensive study of non\-adaptive and residual\-based adaptive sampling for physics\-informed neural networks, Computer Methods in Applied Mechanics and Engineering 403 \(2023\) 115671\.
- \[35\]W\. Gao, C\. Wang, Active learning based sampling for high\-dimensional nonlinear partial differential equations, Journal of Computational Physics 475 \(2023\) 111848\.
- \[36\]G\. K\. R\. Lau, A\. Hemachandra, S\.\-K\. Ng, B\. K\. H\. Low, PINNACLE: PINN adaptive collocation and experimental points selection, in: The Twelfth International Conference on Learning Representations, 2024\.
- \[37\]M\. Tancik, P\. Srinivasan, B\. Mildenhall, S\. Fridovich\-Keil, N\. Raghavan, U\. Singhal, R\. Ramamoorthi, J\. Barron, R\. Ng, Fourier features let networks learn high frequency functions in low dimensional domains, Advances in Neural Information Processing Systems 33 \(2020\) 7537–7547\.
- \[38\]S\. Wang, H\. Wang, P\. Perdikaris, On the eigenvector bias of Fourier feature networks: From regression to solving multi\-scale PDEs with physics\-informed neural networks, Computer Methods in Applied Mechanics and Engineering 384 \(2021\) 113938\.
- \[39\]V\. Sitzmann, J\. Martel, A\. Bergman, D\. Lindell, G\. Wetzstein, Implicit neural representations with periodic activation functions, Advances in neural information processing systems 33 \(2020\) 7462–7473\.
- \[40\]C\. Si, M\. Yan, X\. Li, Z\. Xia,[Complex physics\-informed neural network](https://arxiv.org/abs/2502.04917)\(2025\)\.[arXiv:2502\.04917](http://arxiv.org/abs/2502.04917)\. URL[https://arxiv\.org/abs/2502\.04917](https://arxiv.org/abs/2502.04917)
- \[41\]X\. Zeng, Y\. Zhu, A nurbs\-based parameterization physics\-informed neural network with an adaptive architecture for solving pdes, Journal of Computational Physics 562 \(2026\) 114980\.
- \[42\]Z\. Zhao, X\. Ding, B\. A\. Prakash, PINNsFormer: A transformer\-based framework for physics\-informed neural networks, in: The Twelfth International Conference on Learning Representations, 2024\.
- \[43\]Y\. Wang, J\. Sun, J\. Bai, C\. Anitescu, M\. S\. Eshaghi, X\. Zhuang, T\. Rabczuk, Y\. Liu, Kolmogorov–arnold\-informed neural network: A physics\-informed deep learning framework for solving forward and inverse problems based on kolmogorov–arnold networks, Computer Methods in Applied Mechanics and Engineering 433 \(2025\) 117518\.
- \[44\]A\. D\. Jagtap, G\. E\. Karniadakis, Extended physics\-informed neural networks \(xpinns\): A generalized space\-time domain decomposition based deep learning framework for nonlinear partial differential equations, Communications in Computational Physics 28 \(5\) \(2020\)\.
- \[45\]S\. Wang, B\. Li, Y\. Chen, P\. Perdikaris, Piratenets: Physics\-informed deep learning with residual adaptive networks, Journal of Machine Learning Research 25 \(402\) \(2024\) 1–51\.
- \[46\]A\. H\. Mustajab, H\. Lyu, Z\. Rizvi, F\. Wuttke, Physics\-informed neural networks for high\-frequency and multi\-scale problems using transfer learning, Applied Sciences 14 \(8\) \(2024\) 3204\.
- \[47\]Y\. Wang, C\.\-Y\. Lai, Multi\-stage neural networks: Function approximator of machine precision, Journal of Computational Physics 504 \(2024\) 112865\.
- \[48\]E\. C\. Cyr, M\. A\. Gulian, R\. G\. Patel, M\. Perego, N\. A\. Trask, Robust training and initialization of deep neural networks: An adaptive basis viewpoint, in: Mathematical and Scientific Machine Learning, PMLR, 2020, pp\. 512–536\.
- \[49\]G\.\-B\. Huang, Q\.\-Y\. Zhu, C\.\-K\. Siew, Extreme learning machine: theory and applications, Neurocomputing 70 \(1\-3\) \(2006\) 489–501\.
- \[50\]V\. Dwivedi, B\. Srinivasan, Physics informed extreme learning machine \(pielm\)–a rapid method for the numerical solution of partial differential equations, Neurocomputing 391 \(2020\) 96–118\.
- \[51\]S\. Dong, Z\. Li, Local extreme learning machines and domain decomposition for solving linear and nonlinear partial differential equations, Computer Methods in Applied Mechanics and Engineering 387 \(2021\) 114129\.
- \[52\]Y\. Wu, M\. Aguiar, K\. H\. Johansson, M\. Barreau, Iterative training of physics\-informed neural networks with fourier\-enhanced features, arXiv preprint arXiv:2510\.19399 \(2025\)\.
- \[53\]K\. Levenberg, A method for the solution of certain non\-linear problems in least squares, Quarterly of applied mathematics 2 \(2\) \(1944\) 164–168\.
- \[54\]D\. W\. Marquardt, An algorithm for least\-squares estimation of nonlinear parameters, Journal of the society for Industrial and Applied Mathematics 11 \(2\) \(1963\) 431–441\.
- \[55\]A\. G\. Wilson, Z\. Hu, R\. Salakhutdinov, E\. P\. Xing, Deep kernel learning, in: Artificial Intelligence and Statistics, PMLR, 2016, pp\. 370–378\.
- \[56\]A\. Rahimi, B\. Recht, Random features for large\-scale kernel machines, Advances in Neural Information Processing Systems 20 \(2007\)\.
- \[57\]D\. P\. Kingma, J\. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412\.6980 \(2014\)\.
- \[58\]R\. Coleman, Calculus on normed vector spaces, Choice Reviews Online \(2013\)\.
- \[59\]M\. Penwarden, A\. D\. Jagtap, S\. Zhe, G\. E\. Karniadakis, R\. M\. Kirby, A unified scalable framework for causal sweeping strategies for physics\-informed neural networks \(PINNs\) and their temporal decompositions, Journal of Computational Physics 493 \(2023\) 112464\.[doi:10\.1016/j\.jcp\.2023\.112464](https://doi.org/10.1016/j.jcp.2023.112464)\.
- \[60\]C\. E\. Shannon, Communication in the presence of noise, Proceedings of the IRE 37 \(1\) \(2006\) 10–21\.
- \[61\]Z\. Hao, J\. Yao, C\. Su, H\. Su, Z\. Wang, F\. Lu, Z\. Xia, Y\. Zhang, S\. Liu, L\. Lu, et al\., Pinnacle: A comprehensive benchmark of physics\-informed neural networks for solving pdes, Advances in Neural Information Processing Systems 37 \(2024\) 76721–76774\.