Feature Interaction Modeling for Physics-Informed Neural Networks and Neural Operators

# Feature Interaction Modeling for Physics-Informed Neural Networks and Neural Operators
Source: [https://arxiv.org/html/2607.28762](https://arxiv.org/html/2607.28762)
###### Abstract

This work embeds feature interaction modules derived from factorization machines \(FMs\) into physics\-informed neural networks \(PINNs\) and neural operator learning, to enhance model expressiveness for solution manifolds of parameterized partial differential equations \(PDEs\)\. Motivated by the second\-order Taylor expansion of multivariate functions to characterize variable couplings, we first propose FM\-PINN\. It explicitly captures spatio\-temporal variable interactions and improves the approximation accuracy for smooth high\-order PDEs\. We further group spatial coordinates, time, physical parameters, and initial and boundary conditions into independent feature sets and model their cross\-group interactions\. Based on this strategy, we develop FM\-Operator and FM\-DeepONet, which are particularly effective for nonlinear conservation laws and problems with sharp gradients or discontinuities, while offering no consistent advantage on smooth operator learning benchmarks\. Numerical tests demonstrate that the proposed mechanism delivers substantial accuracy gains on challenging shock\-dominated equations, indicating a promising direction for physics\-consistent modeling of parameterized PDEs with strong cross\-field dependencies\.

## 1Introduction

Partial differential equations \(PDEs\) are the mathematical language of the physical world and play an indispensable role across science and engineering\. Yet solving them efficiently and accurately remains a fundamental challenge\. Traditional numerical methods typically require predefined meshes, spatial discretization, and problem\-specific numerical schemes\. In 2019,Raissiet al\.\([2019](https://arxiv.org/html/2607.28762#bib.bib1)\)introduced physics\-informed neural networks \(PINNs\), a breakthrough that offered a mesh\-free alternative by encoding the governing equations, initial conditions, and boundary conditions directly into the neural network training objective\.

Around the same time,Luet al\.\([2021](https://arxiv.org/html/2607.28762#bib.bib6)\)opened a new frontier by bringing the concept of learning mappings between function spaces into neural PDE solvers\. Their Deep Operator Network \(DeepONet\) was designed to handle families of PDEs—varying initial conditions, boundary conditions, or equation parameters—by representing the solution operator through a branch–trunk factorization\. The branch network encodes an input function sampled at a fixed set of sensors, while the trunk network encodes the query coordinates at which the solution is evaluated\. Building on this operator perspective,Liet al\.\([2021](https://arxiv.org/html/2607.28762#bib.bib7)\)later proposed the Fourier neural operator \(FNO\), which parameterizes integral kernels in the Fourier domain and substantially improved approximation accuracy, though it is most naturally suited to regular grids\. These two architectures—DeepONet and FNO—have since become the dominant paradigms in neural operator learning\.

Input feature design is equally critical for both PINNs and neural operators\. The Fourier feature mapping introduced byTanciket al\.\([2020](https://arxiv.org/html/2607.28762#bib.bib4)\)has been widely adopted to alleviate the spectral bias of multilayer perceptrons\. Building on this,Wanget al\.\([2021](https://arxiv.org/html/2607.28762#bib.bib15)\)applied Fourier features to PINNs and demonstrated improved approximation of high\-frequency and multiscale PDE solutions\. More recently,Fazlianiet al\.\([2025](https://arxiv.org/html/2607.28762#bib.bib5)\)proposed SAFE\-NET, which combines Fourier\-based feature engineering with a single\-hidden\-layer network and an effective optimization strategy, simultaneously improving both the accuracy and training efficiency of PINNs\.

However, existing work largely focuses on domain decomposition\(Jagtap and Karniadakis[2020](https://arxiv.org/html/2607.28762#bib.bib2); Huet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib3); Moseleyet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib18)\), spectral factorization\(Liet al\.[2021](https://arxiv.org/html/2607.28762#bib.bib7); Rahmanet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib19); Tranet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib20); Kossaifiet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib17)\), multimodal fusion\(Zhuet al\.[2025](https://arxiv.org/html/2607.28762#bib.bib13); Liuet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib21)\), and input encoding\(Tanciket al\.[2020](https://arxiv.org/html/2607.28762#bib.bib4); Wanget al\.[2021](https://arxiv.org/html/2607.28762#bib.bib15); Fazlianiet al\.[2025](https://arxiv.org/html/2607.28762#bib.bib5)\), with far less attention paid to interactions among the input features themselves\. Yet the widespread use of CP decomposition\(Choet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib14)\), Tucker decomposition\(Kossaifiet al\.[2023](https://arxiv.org/html/2607.28762#bib.bib17)\), and other tensor decompositions\(Chenet al\.[2025](https://arxiv.org/html/2607.28762#bib.bib23)\)in PINNs and neural operator learning already hints at the importance of cross\-feature interactions\. Motivated by this observation, we introduce factorization machines into PINN and operator learning\. Factorization machines provide an efficient way to model pairwise feature interactions by sharing latent factors across all pairs rather than assigning an independent parameter to each one, thus remaining scalable as the input dimension grows\(Rendle[2010](https://arxiv.org/html/2607.28762#bib.bib9)\)\. Building on this core idea, several variants have since been developed, including NFM\(He and Chua[2017](https://arxiv.org/html/2607.28762#bib.bib10)\), AFM\(Xiaoet al\.[2017](https://arxiv.org/html/2607.28762#bib.bib11)\), and xDeepFM\(Lianet al\.[2018](https://arxiv.org/html/2607.28762#bib.bib12)\)\. In this work, we adapt the bi\-interaction mechanism to continuous PDE inputs and use it as a shared architectural component for both PINNs and neural operators\.

From this formulation, we develop three architectures\.

1. 1\.FM\-PINNexplicitly models pairwise relationships between spatial coordinates and between space and time, outperforming the vanilla PINN on several high\-dimensional smooth problems\.
2. 2\.FM\-Operatoris a new operator architecture that explicitly models the relationships among sensor\-level input functions, query coordinates across multiple dimensions, and auxiliary statistics, achieving superior accuracy on sharp\-gradient problems with fewer parameters \(Figure[1](https://arxiv.org/html/2607.28762#S4.F1)\)\.
3. 3\.FM\-DeepONetinjects field\-interaction blocks into the DeepONet framework to refine its branch and trunk features, substantially improving solution accuracy on problems with sharp gradients and discontinuities\.

## 2Related Work

##### Physics\-informed Neural Networks

To improve the capabilities of PINNs, researchers have addressed their limitations from several directions\.Jagtap and Karniadakis \([2020](https://arxiv.org/html/2607.28762#bib.bib2)\)proposed XPINN, which partitions the space\-time domain into subdomains and assigns an independent network to each, improving local representational capacity and parallelizability, though it typically requires additional interface losses to enforce consistency across subdomains\. Building on this,Huet al\.\([2023](https://arxiv.org/html/2607.28762#bib.bib3)\)introduced APINN, which employs a trainable gating network for soft domain decomposition combined with parameter sharing, thereby partially alleviating the interface\-matching and data\-scarcity issues of hard partitioning\.Choet al\.\([2023](https://arxiv.org/html/2607.28762#bib.bib14)\)proposed SPINN, which uses dimension\-wise separable networks and forward\-mode automatic differentiation to substantially accelerate training and improve accuracy on smooth PDEs\.

##### Neural Operator

For neural operators,Liet al\.\([2023](https://arxiv.org/html/2607.28762#bib.bib16)\)proposed PINO, which incorporates physical constraints into operator learning, combining the generalization of neural operators with the physical consistency of PINNs\. To handle sharp gradients,Lanthaleret al\.\([2023](https://arxiv.org/html/2607.28762#bib.bib8)\)developed Shift\-DeepONet, which uses nonlinear reconstruction to better represent solution operators with discontinuities\. To address the high training cost of FNO,Kossaifiet al\.\([2023](https://arxiv.org/html/2607.28762#bib.bib17)\)proposed MG\-TFNO, which applies Tucker decomposition to the spectral weight tensors of FNO, substantially reducing the parameter count and memory overhead of its Fourier layers\. More recently, research on universal PDE solvers has also gained momentum;Zhuet al\.\([2025](https://arxiv.org/html/2607.28762#bib.bib13)\)proposed PI\-MFM, which converts PDEs into Polish expressions and automatically assembles vector derivative computations from input expressions into PDE residual losses to enable multi\-operator learning\.

## 3Problem Setting

##### Physics\-informed Neural Networks \(PINNs\)

Consider a PDE system:

\{ℱ​\[u\]​\(t,𝐱\)=0,\(t,𝐱\)∈\(0,T\]×Ω,u​\(0,𝐱\)=g​\(𝐱\),𝐱∈Ω,ℬ​\[u\]​\(t,𝐱\)=h​\(t,𝐱\),\(t,𝐱\)∈\(0,T\]×∂Ω\.\\begin\{cases\}\\mathcal\{F\}\[u\]\(t,\\mathbf\{x\}\)=0,&\(t,\\mathbf\{x\}\)\\in\(0,T\]\\times\\Omega,\\\\ u\(0,\\mathbf\{x\}\)=g\(\\mathbf\{x\}\),&\\mathbf\{x\}\\in\\Omega,\\\\ \\mathcal\{B\}\[u\]\(t,\\mathbf\{x\}\)=h\(t,\\mathbf\{x\}\),&\(t,\\mathbf\{x\}\)\\in\(0,T\]\\times\\partial\\Omega\.\\end\{cases\}whereΩ\\Omegadenotes the spatial domain,ℱ​\[⋅\]\\mathcal\{F\}\[\\cdot\]is the PDE residual operator,ℬ​\[⋅\]\\mathcal\{B\}\[\\cdot\]represents the boundary operator, andggandhhdenote the initial and boundary conditions, respectively\.

A PINN approximatesuuwith a neural networkuθ​\(t,𝐱\)u\_\{\\theta\}\(t,\\mathbf\{x\}\)and minimizes

ℒPINN​\(θ\)=λr​ℒr​\(θ\)\+λb​ℒb​\(θ\)\+λ0​ℒ0​\(θ\),\\mathcal\{L\}\_\{\\mathrm\{PINN\}\}\(\\theta\)=\\lambda\_\{r\}\\mathcal\{L\}\_\{r\}\(\\theta\)\+\\lambda\_\{b\}\\mathcal\{L\}\_\{b\}\(\\theta\)\+\\lambda\_\{0\}\\mathcal\{L\}\_\{0\}\(\\theta\),\(1\)whereℒr,ℒb,ℒ0\\mathcal\{L\}\_\{r\},\\mathcal\{L\}\_\{b\},\\mathcal\{L\}\_\{0\}are the mean squared errors of the PDE residual, boundary conditions, and initial conditions, respectively\.

##### Deep Operator Network \(DeepONet\)

Given an input function sampled at sensorsss, an optional physical parameterη\\eta, and a queryyy, operator learning approximates𝒢:\(s,η\)↦u\\mathcal\{G\}:\(s,\\eta\)\\mapsto u\. A standard DeepONet predicts

𝒢θ​\(s,η\)​\(y\)=∑k=1pbk​\(s,η\)​ϕk​\(y\)\+β,\\mathcal\{G\}\_\{\\theta\}\(s,\\eta\)\(y\)=\\sum\_\{k=1\}^\{p\}b\_\{k\}\(s,\\eta\)\\,\\phi\_\{k\}\(y\)\+\\beta,\(2\)and is trained by minimizing the supervised loss

ℒ​\(θ\)=1N​P​∑i=1N∑j=1P\|𝒢θ​\(si,ηi\)​\(yj\)−𝒢​\(si,ηi\)​\(yj\)\|2,\\mathcal\{L\}\(\\theta\)=\\frac\{1\}\{NP\}\\sum\_\{i=1\}^\{N\}\\sum\_\{j=1\}^\{P\}\\bigl\|\\mathcal\{G\}\_\{\\theta\}\(s\_\{i\},\\eta\_\{i\}\)\(y\_\{j\}\)\-\\mathcal\{G\}\(s\_\{i\},\\eta\_\{i\}\)\(y\_\{j\}\)\\bigr\|^\{2\},\(3\)overNNinput functions andPPquery points\.

## 4Continuous\-Field Factorized Interactions

### Feature Interaction Block

Let a PDE\-learning input be partitioned intoMMsemantic fields,

F=\(F1,…,FM\),Fi∈ℝdi\.F=\(F\_\{1\},\\ldots,F\_\{M\}\),\\qquad F\_\{i\}\\in\\mathbb\{R\}^\{d\_\{i\}\}\.Each field has its own affine projection into a sharedrr\-dimensional space,

ei=Wi​Fi\+βi,ei∈ℝr\.e\_\{i\}=W\_\{i\}F\_\{i\}\+\\beta\_\{i\},\\qquad e\_\{i\}\\in\\mathbb\{R\}^\{r\}\.Following the bi\-interaction pooling of neural factorization machines \(NFM\)\(He and Chua[2017](https://arxiv.org/html/2607.28762#bib.bib10)\), we aggregate pairwise field interactions as

zBI\\displaystyle z\_\{\\mathrm\{BI\}\}=∑i<jei⊙ej\\displaystyle=\\sum\_\{i<j\}e\_\{i\}\\odot e\_\{j\}\(4\)=12​\[\(∑iei\)⊙2−∑iei⊙2\],\\displaystyle=\\frac\{1\}\{2\}\\left\[\\left\(\\sum\_\{i\}e\_\{i\}\\right\)^\{\\odot 2\}\-\\sum\_\{i\}e\_\{i\}^\{\\odot 2\}\\right\],where⊙\\odotis the Hadamard product\.

### Proposed Architectures

##### FM\-PINN\.

FM\-PINN treats each coordinate component as a separate field\. For a one\-dimensional time\-dependent PDE,F=\(x,t/T\)F=\(x,\\,t/T\); for a steadyDD\-dimensional problem,F=\(x1,…,xD\)F=\(x\_\{1\},\\dots,x\_\{D\}\)\. Each field is passed through the Feature Interaction Block \(embedding and bi\-interaction pooling as defined in Section 4\.1\), producing the interaction vectorzBIz\_\{\\mathrm\{BI\}\}\. An MLP decoder then predicts the solution:

uθ​\(y\)=MLP​\(zBI\)\.u\_\{\\theta\}\(y\)=\\mathrm\{MLP\}\(z\_\{\\mathrm\{BI\}\}\)\.The network is trained with the same weighted PDE, initial\-condition, and boundary\-condition residuals as the baseline PINN \(Eq\.[1](https://arxiv.org/html/2607.28762#S3.E1)\)\. This definition requiresM≥2M\\geq 2; a one\-field task needs either a different field partition or an explicit first\-order path to avoid a constant prediction\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/FM_Operator.drawio.png)Figure 1:FM\-Operator architecture\.
##### FM\-Operator\.

FM\-Operator is a direct operator architecture that does not rely on a branch–trunk factorization\. The input fields areF=\(s,y,q​\(s\),η\)F=\(s,y,q\(s\),\\eta\), whereq​\(s\)=\[mean⁡\(s\),std⁡\(s\),max⁡\|s\|\]q\(s\)=\[\\operatorname\{mean\}\(s\),\\operatorname\{std\}\(s\),\\max\|s\|\]andη\\etais omitted when already part ofss\. After the Feature Interaction Block produceszBIz\_\{\\mathrm\{BI\}\}, the prediction combines an MLP path with a linear first\-order path:

𝒢θ​\(s,η\)​\(y\)=MLP​\(zBI\)\+Linear​\(\[F1;…;FM\]\)\.\\mathcal\{G\}\_\{\\theta\}\(s,\\eta\)\(y\)=\\mathrm\{MLP\}\(z\_\{\\mathrm\{BI\}\}\)\+\\mathrm\{Linear\}\\bigl\(\[F\_\{1\};\\ldots;F\_\{M\}\]\\bigr\)\.The linear path retains raw first\-order effects, while the MLP path operates on pooled cross\-field products\. FM\-Operator is evaluated in two granularity variants:*split*, where each query coordinate is its own field \(y1,…,ydyy\_\{1\},\\dots,y\_\{d\_\{y\}\}\), and*non\-split*, where the full query vectoryyis a single field\.

##### FM\-DeepONet\.

FM\-DeepONet injects the Feature Interaction Block into the standard DeepONet framework\. Letb,ϕ∈ℝpb,\\phi\\in\\mathbb\{R\}^\{p\}denote the branch and trunk features produced by the baseline DeepONet encoders\. We construct the field setF=\(b,ϕ,q​\(s\),η\)F=\(b,\\phi,q\(s\),\\eta\)and apply the Feature Interaction Block to obtainzBIz\_\{\\mathrm\{BI\}\}\. An MLP decodeszBIz\_\{\\mathrm\{BI\}\}into a refined representation, which is then projected onto correction terms:

Δ​b=Pb​\(r\),Δ​ϕ=Pϕ​\(r\),\\Delta b=P\_\{b\}\(r\),\\qquad\\Delta\\phi=P\_\{\\phi\}\(r\),wherer=MLP​\(zBI\)r=\\mathrm\{MLP\}\(z\_\{\\mathrm\{BI\}\}\)andPb,PϕP\_\{b\},P\_\{\\phi\}are trainable affine projections\. The corrected featuresb~=b\+Δ​b,ϕ~=ϕ\+Δ​ϕ\\tilde\{b\}=b\+\\Delta b,\\;\\tilde\{\\phi\}=\\phi\+\\Delta\\phireplace the original ones in the DeepONet inner\-product head:

𝒢θ​\(s\)​\(y\)=1p​⟨b~,ϕ~⟩\+β\.\\mathcal\{G\}\_\{\\theta\}\(s\)\(y\)=\\frac\{1\}\{\\sqrt\{p\}\}\\langle\\tilde\{b\},\\tilde\{\\phi\}\\rangle\+\\beta\.This design retains the algebraic structure of DeepONet while allowing pairwise interactions among branch features, trunk features, and auxiliary inputs to jointly refine the prediction\.

### Theoretical Motivation

We motivate the bi\-interaction design by connecting it to a block\-wise second\-order Taylor expansion\. For a scalar targetg​\(F\)g\(F\)that isC3C^\{3\}nearF0F^\{0\}, expanding by field blocks gives

g​\(F0\+δ​F\)=\\displaystyle g\(F^\{0\}\+\\delta F\)=\{\}g​\(F0\)\+∑i∇ig​\(F0\)⊤​δ​Fi\\displaystyle g\(F^\{0\}\)\+\\sum\_\{i\}\\nabla\_\{i\}g\(F^\{0\}\)^\{\\top\}\\delta F\_\{i\}\+12​∑iδ​Fi⊤​Hi​i​\(F0\)​δ​Fi\\displaystyle\+\\frac\{1\}\{2\}\\sum\_\{i\}\\delta F\_\{i\}^\{\\top\}H\_\{ii\}\(F^\{0\}\)\\delta F\_\{i\}\+∑i<jδ​Fi⊤​Hi​j​\(F0\)​δ​Fj\+R3​\(δ​F\),\\displaystyle\+\\sum\_\{i<j\}\\delta F\_\{i\}^\{\\top\}H\_\{ij\}\(F^\{0\}\)\\delta F\_\{j\}\+R\_\{3\}\(\\delta F\),\(5\)whereHi​j=∇Fi​Fj2gH\_\{ij\}=\\nabla^\{2\}\_\{F\_\{i\}F\_\{j\}\}gand the remainder satisfies\|R3\|≤C3​‖δ​F‖3/6\|R\_\{3\}\|\\leq C\_\{3\}\\\|\\delta F\\\|^\{3\}/6when the third derivative is bounded byC3C\_\{3\}\. Equation \([4](https://arxiv.org/html/2607.28762#S4.E4)\) directly exposes the mixed\-block products∑i<j\\sum\_\{i<j\}, matching the cross\-field Hessian terms in \([5](https://arxiv.org/html/2607.28762#S4.E5)\)\. The affine embeddings can represent local centering sinceWi​\(Fi−Fi0\)=Wi​Fi−Wi​Fi0W\_\{i\}\(F\_\{i\}\-F\_\{i\}^\{0\}\)=W\_\{i\}F\_\{i\}\-W\_\{i\}F\_\{i\}^\{0\}\. For aDD\-dimensional coordinate input there are\(D2\)\\binom\{D\}\{2\}such coordinate pairs; for a 1D space–time input the only mixed term is locally proportional toδ​x​δ​t\\delta x\\,\\delta t\.

##### Taylor Expansion for FM\-PINN

For a fixed PDE instance, FM\-PINN learns the coordinate\-to\-solution mapy↦u​\(y\)y\\mapsto u\(y\)\. For a steadyDD\-dimensional problem, each coordinate is treated as one scalar field,Fi=xiF\_\{i\}=x\_\{i\}\. Withδ​xi=xi−xi0\\delta x\_\{i\}=x\_\{i\}\-x\_\{i\}^\{0\}, Equation \([5](https://arxiv.org/html/2607.28762#S4.E5)\) becomes

u​\(x0\+δ​x\)=\\displaystyle u\(x^\{0\}\+\\delta x\)=\{\}u​\(x0\)\+∑i=1Duxi​\(x0\)​δ​xi\\displaystyle u\(x^\{0\}\)\+\\sum\_\{i=1\}^\{D\}u\_\{x\_\{i\}\}\(x^\{0\}\)\\delta x\_\{i\}\+12​∑i=1Duxi​xi​\(x0\)​\(δ​xi\)2\\displaystyle\+\\frac\{1\}\{2\}\\sum\_\{i=1\}^\{D\}u\_\{x\_\{i\}x\_\{i\}\}\(x^\{0\}\)\(\\delta x\_\{i\}\)^\{2\}\+∑i<juxi​xj​\(x0\)​δ​xi​δ​xj\+R3\.\\displaystyle\+\\sum\_\{i<j\}u\_\{x\_\{i\}x\_\{j\}\}\(x^\{0\}\)\\delta x\_\{i\}\\delta x\_\{j\}\+R\_\{3\}\.\(6\)The last sum contains the cross\-coordinate terms explicitly targeted by FM\-PINN\. ADD\-dimensional input has\(D2\)\\binom\{D\}\{2\}such pairs, and the interaction block is most useful when their mixed derivatives are substantial\. For a one\-dimensional time\-dependent problem, the two fields arexxandtt, and the mixed term isux​t​δ​x​δ​tu\_\{xt\}\\delta x\\delta t\. Thus, 1D Burgers still contains a genuine space–time interaction\.

##### Taylor Expansion for FM\-Operator

After sensor discretization, write the target operator as

g​\(s,y,η\)=𝒢​\(s,η\)​\(y\),g\(s,y,\\eta\)=\\mathcal\{G\}\(s,\\eta\)\(y\),\(7\)wheressis the sensor vector,yythe query coordinate, andη\\etaan optional physical parameter\. Around\(s0,y0,η0\)\(s^\{0\},y^\{0\},\\eta^\{0\}\), the second\-order expansion is

g=\\displaystyle g=\{\}g0\+gs⊤​δ​s\+gy⊤​δ​y\+gη⊤​δ​η\\displaystyle g\_\{0\}\+g\_\{s\}^\{\\top\}\\delta s\+g\_\{y\}^\{\\top\}\\delta y\+g\_\{\\eta\}^\{\\top\}\\delta\\eta\+12​δ​s⊤​Hs​s​δ​s\+12​δ​y⊤​Hy​y​δ​y\+12​δ​η⊤​Hη​η​δ​η\\displaystyle\+\\frac\{1\}\{2\}\\delta s^\{\\top\}H\_\{ss\}\\delta s\+\\frac\{1\}\{2\}\\delta y^\{\\top\}H\_\{yy\}\\delta y\+\\frac\{1\}\{2\}\\delta\\eta^\{\\top\}H\_\{\\eta\\eta\}\\delta\\eta\+δ​s⊤​Hs​y​δ​y\+δ​s⊤​Hs​η​δ​η\+δ​y⊤​Hy​η​δ​η\+R3\.\\displaystyle\+\\delta s^\{\\top\}H\_\{sy\}\\delta y\+\\delta s^\{\\top\}H\_\{s\\eta\}\\delta\\eta\+\\delta y^\{\\top\}H\_\{y\\eta\}\\delta\\eta\+R\_\{3\}\.\(8\)FM\-Operator treatsssandyyas separate fields, so its main explicit cross\-field term corresponds toδ​s⊤​Hs​y​δ​y\\delta s^\{\\top\}H\_\{sy\}\\delta y\. Ifη\\etais separate, the two parameter\-related terms are also available; if it is concatenated withss, they belong to the sensor block\.

Linearity inssdoes not remove the sensor–query interaction\. For example, if

𝒢​\(s\)​\(y\)=∑m=1nssm​ϕm​\(y\),\\mathcal\{G\}\(s\)\(y\)=\\sum\_\{m=1\}^\{n\_\{s\}\}s\_\{m\}\\phi\_\{m\}\(y\),\(9\)then∇s​s2g=0\\nabla\_\{ss\}^\{2\}g=0, but

∂2g∂sm​∂yk=∂ϕm​\(y\)∂yk\\frac\{\\partial^\{2\}g\}\{\\partial s\_\{m\}\\partial y\_\{k\}\}=\\frac\{\\partial\\phi\_\{m\}\(y\)\}\{\\partial y\_\{k\}\}\(10\)is generally nonzero\. The relevant question is therefore whether the joint map\(s,y\)↦𝒢​\(s\)​\(y\)\(s,y\)\\mapsto\\mathcal\{G\}\(s\)\(y\)depends strongly on both fields, rather than on whether the underlying PDE is linear\.

##### Inductive bias, not a universal approximation claim\.

A sufficiently large concatenation MLP can approximate multiplication and therefore represent the same cross\-field products\. The proposed block does not prove superior expressiveness; it makes these products explicit before deep nonlinear processing\. Moreover, because the bi\-interaction pooling first projects each field through a sharedrr\-dimensional embedding before computing pairwise Hadamard products, the resulting cross\-field bilinear forms are collectively constrained to rank at mostrr\. Consequently, a collection of arbitrary rank\-rrcross\-field matrices need not be jointly representable\. This reduced parameterization is a compact inductive bias, but it can also cause negative transfer, consistent with the failure cases reported below\.

## 5Experiments

We first describe the task\-specific training setups, then compare FM\-PINN against the vanilla PINN, and finally evaluate the neural operator models\.

### Experimental Setup

##### PINN task\.

PINN applies a linear projection to the raw coordinates, whereas FM\-PINN first maps each coordinate to a 64\-dimensional field embedding and pools its pairwise interactions\. They share the same decoder: a linear64→12864\\\!\\to\\\!128layer, four width\-128 hidden layers with SiLU activations, and a scalar output\. Both models are trained for 30,000 AdamW steps with learning rate10−310^\{\-3\}and weight decay10−610^\{\-6\}\. Latin hypercube sampling givesNf=10,000N\_\{f\}=10\{,\}000interior points andNbc/ic=400N\_\{\\mathrm\{bc/ic\}\}=400boundary or initial points\. Each update uses 1,024 interior points and all 400 boundary/initial points\. The two models use the same samples, residual loss, and evaluation points; Table[2](https://arxiv.org/html/2607.28762#S5.T2)reports one run per model\.

For the steady convection–diffusion, reaction–diffusion, and linear advection tests,Ω=\[0,1\]d\\Omega=\[0,1\]^\{d\}, with the dimensions listed in Table[2](https://arxiv.org/html/2607.28762#S5.T2)\. Their shared manufactured solution is

u∗​\(x\)=∑m=15cm​∏j=1dsin⁡\(km,j​π​xj\),u^\{\*\}\(x\)=\\sum\_\{m=1\}^\{5\}c\_\{m\}\\prod\_\{j=1\}^\{d\}\\sin\(k\_\{m,j\}\\pi x\_\{j\}\),\(11\)wherec=\(1,0\.35,−0\.25,0\.20,0\.10\)c=\(1,0\.35,\-0\.25,0\.20,0\.10\)\. The frequency vectors, truncated to dimensiondd, arek1=\(1,1,…\)k\_\{1\}=\(1,1,\\ldots\),k2=\(2,1,…\)k\_\{2\}=\(2,1,\\ldots\),k3=\(1,2,1,…\)k\_\{3\}=\(1,2,1,\\ldots\),k4=\(2,1,2,1,…\)k\_\{4\}=\(2,1,2,1,\\ldots\), andk5=\(3,1,2,1,…\)k\_\{5\}=\(3,1,2,1,\\ldots\)\. We solve

−0\.05​Δ​u\+β⋅∇u\\displaystyle\-05\\Delta u\+\\beta\\\!\\cdot\\\!\\nabla u=fcd,\\displaystyle=f\_\{\\rm cd\},\(12\)−0\.05​Δ​u\+\(u3−u\)\\displaystyle\-05\\Delta u\+\(u^\{3\}\-u\)=frd,\\displaystyle=f\_\{\\rm rd\},β⋅∇u\\displaystyle\\beta\\\!\\cdot\\\!\\nabla u=fadv,\\displaystyle=f\_\{\\rm adv\},whereβ=\(1,−0\.7,0\.5,−0\.3,0\.2\)\\beta=\(1,\-0\.7,0\.5,\-0\.3,0\.2\)is truncated todd\. Each forcing term is obtained by applying its left\-hand operator to Equation \([11](https://arxiv.org/html/2607.28762#S5.E11)\), and the Dirichlet value isu=u∗u=u^\{\*\}on∂Ω\\partial\\Omega\. Thus the reference values for these tests are evaluated directly fromu∗u^\{\*\}\.

The 4D heat problem is posed on\[0,1\]3×\[0,0\.25\]\[0,1\]^\{3\}\\times\[0,0\.25\]:

ut−0\.05​Δ​u=0,u∗=∑mcm​∏j=13sin⁡\(km,j​π​xj\)×e−0\.05​π2​∥km∥22​t\.\\begin\{gathered\}u\_\{t\}\-0\.05\\Delta u=0,\\\\ u^\{\*\}=\\sum\_\{m\}c\_\{m\}\\prod\_\{j=1\}^\{3\}\\sin\(k\_\{m,j\}\\pi x\_\{j\}\)\\\\ \\times e^\{\-0\.05\\pi^\{2\}\\lVert k\_\{m\}\\rVert\_\{2\}^\{2\}t\}\.\\end\{gathered\}\(13\)Its initial and boundary data are taken from this exact solution\. The 1D Burgers test instead solves

ut\+u​ux−0\.01​ux​x=0,\(x,t\)∈\[−1,1\]×\[0,1\],u\_\{t\}\+uu\_\{x\}\-0\.01u\_\{xx\}=0,\\quad\(x,t\)\\in\[\-1,1\]\\times\[0,1\],\(14\)withu​\(x,0\)=−sin⁡\(π​x\)u\(x,0\)=\-\\sin\(\\pi x\)andu​\(−1,t\)=u​\(1,t\)=0u\(\-1,t\)=u\(1,t\)=0\.

The two remaining 4D time\-dependent tests use\(x,y,z\)∈\[0,1\]3\(x,y,z\)\\in\[0,1\]^\{3\}andt∈\[0,1\]t\\in\[0,1\]\. The Burgers\-type equation is

ut\+u​\(ux\+uy\+uz\)−0\.01​Δ​u=f,u\_\{t\}\+u\(u\_\{x\}\+u\_\{y\}\+u\_\{z\}\)\-0\.01\\Delta u=f,\(15\)with manufactured solution

u∗=∑mcm​∏j=13sin⁡\(km,j​π​xj\)​e−0\.01​π2​∥km∥22​t,u^\{\*\}=\\sum\_\{m\}c\_\{m\}\\prod\_\{j=1\}^\{3\}\\sin\(k\_\{m,j\}\\pi x\_\{j\}\)e^\{\-0\.01\\pi^\{2\}\\lVert k\_\{m\}\\rVert\_\{2\}^\{2\}t\},\(16\)where\(km,cm\)\(k\_\{m\},c\_\{m\}\)are\(\(1,1,1\),−0\.90\)\(\(1,1,1\),\-0\.90\),\(\(2,1,1\),0\.25\)\(\(2,1,1\),0\.25\),\(\(1,2,1\),−0\.20\)\(\(1,2,1\),\-0\.20\), and\(\(1,1,2\),0\.15\)\(\(1,1,2\),0\.15\)\. We setf=ut∗\+u∗​\(ux∗\+uy∗\+uz∗\)−0\.01​Δ​u∗f=u\_\{t\}^\{\*\}\+u^\{\*\}\(u\_\{x\}^\{\*\}\+u\_\{y\}^\{\*\}\+u\_\{z\}^\{\*\}\)\-0\.01\\Delta u^\{\*\}, and impose the initial and boundary data fromu∗u^\{\*\}\. The wave problem is

ut​t−Δ​u=0,u∗=∑mcm​∏j=13sin⁡\(km,j​π​xj\)×cos⁡\(π​∥km∥2​t\),\\begin\{gathered\}u\_\{tt\}\-\\Delta u=0,\\\\ u^\{\*\}=\\sum\_\{m\}c\_\{m\}\\prod\_\{j=1\}^\{3\}\\sin\(k\_\{m,j\}\\pi x\_\{j\}\)\\\\ \\times\\cos\(\\pi\\lVert k\_\{m\}\\rVert\_\{2\}t\),\\end\{gathered\}\(17\)The five\(km,cm\)\(k\_\{m\},c\_\{m\}\)pairs are\(\(1,1,1\),1\)\(\(1,1,1\),1\),\(\(2,1,1\),0\.35\)\(\(2,1,1\),0\.35\),\(\(1,2,3\),−0\.25\)\(\(1,2,3\),\-0\.25\),\(\(2,2,1\),0\.20\)\(\(2,2,1\),0\.20\), and\(\(3,1,2\),0\.10\)\(\(3,1,2\),0\.10\)\. Its initial displacement and boundary values come fromu∗u^\{\*\}, andut​\(x,y,z,0\)=0u\_\{t\}\(x,y,z,0\)=0\.

##### Neural operator task\.

The compared models are listed in Table[1](https://arxiv.org/html/2607.28762#S5.T1)\. In the split/non\-split study and the four\-equation comparison, each nonlinear dataset contains 256 training, 64 validation, and 64 test functions\. An input function is sampled at 64 sensors, and each solution is stored on a128×80128\\times 80space–time grid\. Every AdamW update samples 32 functions and 128 query points per function\. Training lasts 30,000 steps with learning rate10−310^\{\-3\}, weight decay10−610^\{\-6\}, and SiLU activations\. DeepONet and FM\-DeepONet use width\-128, depth\-3 branch and trunk networks with a 128\-dimensional latent basis\. The NFM embedding and interaction decoder also have width 128\. We report the mean and sample standard deviation of the global relativeL2L\_\{2\}error\. \(In the split/non\-split study, the reference solutions are generated by a periodic finite\-volume solver with Rusanov flux and then downsampled to the learning grid\.\)

The cubic conservation and Buckley–Leverett equations are

ut\+\(u3/3\)x=0\.002​ux​x,ut\+\[u2u2\+0\.5​\(1−u\)2\]x=0,\\begin\{gathered\}u\_\{t\}\+\(u^\{3\}/3\)\_\{x\}=0\.002u\_\{xx\},\\\\ u\_\{t\}\+\\left\[\\frac\{u^\{2\}\}\{u^\{2\}\+0\.5\(1\-u\)^\{2\}\}\\right\]\_\{x\}=0,\\end\{gathered\}\(18\)on the periodic intervalx∈\[0,1\)x\\in\[0,1\), witht∈\[0,1\]t\\in\[0,1\]andt∈\[0,0\.6\]t\\in\[0,0\.6\], respectively\. Their six\-mode random Fourier initial fields have coefficient magnitudes decaying as1/k1/k\. Cubic fields have zero offset and amplitude 0\.9; Buckley–Leverett fields have offset 0\.5 and amplitude 0\.42 and are clipped to\[0\.02,0\.98\]\[0\.02,0\.98\]\. For the four\-equation multi\-model comparison, reference solutions are computed by WENO5 with SSPRK3 \(512 spatial points for cubic conservation; 1,024 for Buckley–Leverett\) and downsampled to the learning grid\.

Kuramoto–Sivashinsky solves

ut\+u​ux\+ux​x\+ux​x​x​x=0u\_\{t\}\+uu\_\{x\}\+u\_\{xx\}\+u\_\{xxxx\}=0\(19\)on a periodic interval of length 22 fort∈\[0,1\]t\\in\[0,1\]\. Its initial fields use the same six\-mode1/k1/kconstruction with amplitude 0\.6\. The reference is generated by a Fourier pseudospectral ETDRK4 solver with time step 0\.002 and two\-thirds dealiasing\. For LWR traffic flow,

ρt\+\[ρ​\(1−ρ\)\]x=0,\(x,t\)∈\[0,1\]×\[0,0\.8\]\.\\rho\_\{t\}\+\[\\rho\(1\-\\rho\)\]\_\{x\}=0,\\qquad\(x,t\)\\in\[0,1\]\\times\[0,0\.8\]\.\(20\)Each input is a Riemann state withρL∼𝒰​\(0\.05,0\.45\)\\rho\_\{L\}\\sim\\mathcal\{U\}\(0\.05,0\.45\),ρR∼𝒰​\(0\.55,0\.95\)\\rho\_\{R\}\\sim\\mathcal\{U\}\(0\.55,0\.95\), and initial jumpx0∼𝒰​\(0\.30,0\.70\)x\_\{0\}\\sim\\mathcal\{U\}\(0\.30,0\.70\)\. Its analytic reference places the jump atx0\+\(1−ρL−ρR\)​tx\_\{0\}\+\(1\-\\rho\_\{L\}\-\\rho\_\{R\}\)t\.

Square\-pulse advection solvesut\+0\.6​ux=0u\_\{t\}\+0\.6u\_\{x\}=0on the periodic domain\[0,1\)×\[0,1\]\[0,1\)\\times\[0,1\]\. The pulse amplitude, width, and center followA∼𝒰​\(0\.6,1\.4\)A\\sim\\mathcal\{U\}\(0\.6,1\.4\),w∼𝒰​\(0\.12,0\.28\)w\\sim\\mathcal\{U\}\(0\.12,0\.28\), andxc∼𝒰​\(0\.20,0\.80\)x\_\{c\}\\sim\\mathcal\{U\}\(0\.20,0\.80\)\. Its exact reference is the periodic translationxc​\(t\)=\(xc\+0\.6​t\)mod1x\_\{c\}\(t\)=\(x\_\{c\}\+0\.6t\)\\bmod 1\.

The periodic Burgers benchmark solves

ut\+u​ux=0\.01​ux​x,\(x,t\)∈\[0,1\)×\[0,1\],u\_\{t\}\+uu\_\{x\}=0\.01u\_\{xx\},\\qquad\(x,t\)\\in\[0,1\)\\times\[0,1\],\(21\)from normalized six\-mode periodic Fourier fields\. Its reference uses a periodic finite\-volume solver with Rusanov flux, centered diffusion, and time step5×10−45\\times 10^\{\-4\}\. In the parametric Burgers benchmark, the same equation is posed on\[−1,1\]×\[0,1\]\[\-1,1\]\\times\[0,1\], withν∼𝒰​\(0\.005,0\.02\)\\nu\\sim\\mathcal\{U\}\(0\.005,0\.02\)and Dirichlet valuesbL,bR∼𝒰​\(−0\.3,0\.3\)b\_\{L\},b\_\{R\}\\sim\\mathcal\{U\}\(\-0\.3,0\.3\)\. Withz=\(x\+1\)/2z=\(x\+1\)/2,

u​\(x,0\)=\(1−z\)​bL\+z​bR\+∑k=16ak​sin⁡\(k​π​z\),u\(x,0\)=\(1\-z\)b\_\{L\}\+zb\_\{R\}\+\\sum\_\{k=1\}^\{6\}a\_\{k\}\\sin\(k\\pi z\),\(22\)whereaka\_\{k\}decays as1/k1/kand the perturbation is scaled to amplitude at most 0\.8\. The numerical reference uses upwind convection, centered diffusion, and time step5×10−45\\times 10^\{\-4\}, enforcingbL,bRb\_\{L\},b\_\{R\}at every step\. The branch input is the 64 initial\-condition sensors followed by normalizedν,bL,bR\\nu,b\_\{L\},b\_\{R\}, for 67 entries\.

Finally, the separate Shift\-DeepONet comparison uses 256 training and 64 test functions, 64 input sensors, and a512×80512\\times 80output grid\. It solvesut\+u​ux=0\.001​ux​xu\_\{t\}\+uu\_\{x\}=0\.001u\_\{xx\}on\[−1,1\]×\[0,1\]\[\-1,1\]\\times\[0,1\], with homogeneous Dirichlet boundaries and random six\-mode sine initial data\. The reference uses conservative WENO5, SSPRK3, and maximum time step2×10−42\\times 10^\{\-4\}\. FM\-Operator and Shift\-DeepONet share the same three seeds, functions, batches, optimization budget, and checkpoint rule\. Compared to earlier versions, the reference solutions for all operator benchmarks have been redesigned with higher\-resolution numerical schemes to ensure reliable error evaluation\.

ModelParametersDeepONet∼\\sim108KDeepONet \(short\)∼\\sim58KFM\-DeepONet \(short\)∼\\sim125KFM\-DeepONet∼\\sim208KShift\-DeepONet∼\\sim141KFM\-Operator∼\\sim43KTable 1:Parameter counts for the compared neural operator models\.PDEPINNFM\-PINNError reductionPoisson \(2D\)1\.850×10−21\.850\{\\times\}10^\{\-2\}1\.782×10−21\.782\{\\times\}10^\{\-2\}3\.64%3\.64\\%Darcy \(2D\)2\.411×10−32\.411\{\\times\}10^\{\-3\}3\.156×10−33\.156\{\\times\}10^\{\-3\}−30\.91%\-30\.91\\%Burgers \(1D,ν=0\.01\\nu=0\.01\)1\.249×10−21\.249\{\\times\}10^\{\-2\}1\.237×10−21\.237\{\\times\}10^\{\-2\}0\.96%0\.96\\%Heat \(4D\)8\.018×10−38\.018\{\\times\}10^\{\-3\}5\.014×10−35\.014\{\\times\}10^\{\-3\}37\.46%37\.46\\%Poisson \(4D\)1\.352×10−11\.352\{\\times\}10^\{\-1\}1\.298×10−11\.298\{\\times\}10^\{\-1\}3\.99%3\.99\\%Poisson \(10D\)4\.589×10−14\.589\{\\times\}10^\{\-1\}4\.225×10−14\.225\{\\times\}10^\{\-1\}7\.92%7\.92\\%Convection–Diffusion \(2D\)1\.898×10−31\.898\{\\times\}10^\{\-3\}1\.031×10−31\.031\{\\times\}10^\{\-3\}45\.69%45\.69\\%Convection–Diffusion \(4D\)2\.602×10−22\.602\{\\times\}10^\{\-2\}1\.458×10−21\.458\{\\times\}10^\{\-2\}43\.98%43\.98\\%Convection–Diffusion \(5D\)3\.817×10−23\.817\{\\times\}10^\{\-2\}3\.012×10−23\.012\{\\times\}10^\{\-2\}21\.08%21\.08\\%Reaction–Diffusion \(2D\)4\.608×10−44\.608\{\\times\}10^\{\-4\}5\.672×10−45\.672\{\\times\}10^\{\-4\}−23\.09%\-23\.09\\%Reaction–Diffusion \(4D\)1\.498×10−21\.498\{\\times\}10^\{\-2\}1\.272×10−21\.272\{\\times\}10^\{\-2\}15\.14%15\.14\\%Reaction–Diffusion \(5D\)3\.270×10−23\.270\{\\times\}10^\{\-2\}2\.632×10−22\.632\{\\times\}10^\{\-2\}19\.51%19\.51\\%Linear Advection \(2D\)4\.213×10−34\.213\{\\times\}10^\{\-3\}2\.495×10−32\.495\{\\times\}10^\{\-3\}40\.78%40\.78\\%Linear Advection \(4D\)2\.111×10−22\.111\{\\times\}10^\{\-2\}7\.755×10−37\.755\{\\times\}10^\{\-3\}63\.26%63\.26\\%Linear Advection \(5D\)3\.133×10−23\.133\{\\times\}10^\{\-2\}1\.396×10−21\.396\{\\times\}10^\{\-2\}55\.43%55\.43\\%Darcy \(4D\)1\.066×10−11\.066\{\\times\}10^\{\-1\}1\.005×10−11\.005\{\\times\}10^\{\-1\}5\.73%5\.73\\%Wave \(4D\)4\.988×10−14\.988\{\\times\}10^\{\-1\}3\.149×10−13\.149\{\\times\}10^\{\-1\}36\.87%36\.87\\%Burgers \(4D, smooth\)6\.489×10−36\.489\{\\times\}10^\{\-3\}3\.886×10−33\.886\{\\times\}10^\{\-3\}40\.12%40\.12\\%Table 2:RelativeL2L^\{2\}error comparison between PINN and FM\-PINN\. The error reduction is computed as\(PINN−FM​\-​PINN\)/PINN×100%\(\\mathrm\{PINN\}\-\\mathrm\{FM\}\\mbox\{\-\}\\mathrm\{PINN\}\)/\\mathrm\{PINN\}\\times 100\\%; negative values indicate that FM\-PINN has a higher error than PINN\.

### Result Analysis on FM\-PINN vs\. PINN

Table[2](https://arxiv.org/html/2607.28762#S5.T2)mainly reports equations with two or more spatial coordinates, where interactions between coordinates are significant\. We also keep 1D Burgers because its input includes both space and time,\(x,t\)\(x,t\), and the solution develops a steep nonlinear profile\.

FM\-PINN achieves a lower error on 16 of the 18 cases\. The improvement is more pronounced in several higher\-dimensional problems\. For example, as depicted in Figure[2](https://arxiv.org/html/2607.28762#S5.F2), the error is reduced by 40\.78% for 2D linear advection, 43\.98% for 4D convection–diffusion, 55\.43% for 5D linear advection, and 40\.12% for 4D Burgers\. The two exceptions are 2D Darcy and 2D reaction–diffusion, where the error increases by 30\.91% and 23\.09%, respectively\. Overall, the results show that feature interactions are most useful as the number of input coordinates grows\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/linear_advection2d_mlp_pinn_best_triptych.png)

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/linear_advection2d_nfm_pinn_best_triptych.png)

Figure 2:Comparison of PINN \(top\) and FM\-PINN \(bottom\) on 2D linear advection\.
### Result Analysis on Neural Operators

##### FM\-Operator \(split\) vs\. FM\-Operator \(non\-split\)\.

The split model treats each query coordinate as a separate field, while the non\-split model keeps the full query vector in one field\. Their errors are close on most benchmarks\. The difference is only 0\.05% on parametric Burgers, 1\.21% on Buckley–Leverett, 1\.52% on LWR, and 2\.59% on square\-pulse advection\. Splitting the query is therefore not the main source of the gain, as depicted in Table[3](https://arxiv.org/html/2607.28762#S5.T3); interactions between coordinates inside the query contribute less than interactions between the query and the input function\. Compared with the approximately 108K\-parameter DeepONet, the best FM\-Operator variant reduces the mean error by 17\.28% on Buckley–Leverett, 54\.37% on parametric Burgers \(Figure[3](https://arxiv.org/html/2607.28762#S5.F3)\), 38\.06% on cubic conservation, 28\.71% on LWR, 36\.96% on periodic Burgers, and 13\.51% on square\-pulse advection\. DeepONet remains stronger on many smooth problems, especially Helmholtz, Poisson, and Heat; the full smooth\-problem results are given in the supplementary material\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/burgers1d_parametric_bcnu_deeponet_best_trainloss_test0.png)

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/burgers1d_parametric_bcnu_nfm_deeponet_best_trainloss_test0.png)

Figure 3:Comparison of DeepONet \(top\) and FM\-Operator \(bottom\) on parametric Burgers\.EquationFM\-Operator\(NFM\-Split\)FM\-Operator\(NFM\-Non\-Split\)DeepONetBuckley–Leverett 1D\(7\.531×10−2\)±\(2\.981×10−3\)\\color\[rgb\]\{0,0,1\}\{\(7\.531\\times 10^\{\-2\}\)\\pm\(2\.981\\times 10^\{\-3\}\)\}\(7\.622×10−2\)±\(5\.044×10−3\)\(7\.622\\times 10^\{\-2\}\)\\pm\(5\.044\\times 10^\{\-3\}\)\(9\.104×10−2\)±\(2\.861×10−3\)\(9\.104\\times 10^\{\-2\}\)\\pm\(2\.861\\times 10^\{\-3\}\)Parametric Burgers 1D\(3\.896×10−2\)±\(5\.403×10−3\)\(3\.896\\times 10^\{\-2\}\)\\pm\(5\.403\\times 10^\{\-3\}\)\(3\.894×10−2\)±\(6\.682×10−3\)\\color\[rgb\]\{0,0,1\}\{\(3\.894\\times 10^\{\-2\}\)\\pm\(6\.682\\times 10^\{\-3\}\)\}\(8\.534×10−2\)±\(8\.611×10−3\)\(8\.534\\times 10^\{\-2\}\)\\pm\(8\.611\\times 10^\{\-3\}\)Cubic conservation 1D\(1\.450×10−1\)±\(1\.944×10−2\)\\color\[rgb\]\{0,0,1\}\{\(1\.450\\times 10^\{\-1\}\)\\pm\(1\.944\\times 10^\{\-2\}\)\}\(1\.559×10−1\)±\(1\.302×10−2\)\(1\.559\\times 10^\{\-1\}\)\\pm\(1\.302\\times 10^\{\-2\}\)\(2\.341×10−1\)±\(1\.616×10−2\)\(2\.341\\times 10^\{\-1\}\)\\pm\(1\.616\\times 10^\{\-2\}\)Kuramoto–Sivashinsky 1D\(2\.834×10−2\)±\(1\.816×10−3\)\(2\.834\\times 10^\{\-2\}\)\\pm\(1\.816\\times 10^\{\-3\}\)\(2\.762×10−2\)±\(3\.950×10−3\)\(2\.762\\times 10^\{\-2\}\)\\pm\(3\.950\\times 10^\{\-3\}\)\(1\.830×10−2\)±\(1\.627×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.830\\times 10^\{\-2\}\)\\pm\(1\.627\\times 10^\{\-3\}\)\}LWR 1D\(5\.793×10−2\)±\(5\.200×10−3\)\\color\[rgb\]\{0,0,1\}\{\(5\.793\\times 10^\{\-2\}\)\\pm\(5\.200\\times 10^\{\-3\}\)\}\(5\.881×10−2\)±\(4\.556×10−3\)\(5\.881\\times 10^\{\-2\}\)\\pm\(4\.556\\times 10^\{\-3\}\)\(8\.126×10−2\)±\(3\.710×10−3\)\(8\.126\\times 10^\{\-2\}\)\\pm\(3\.710\\times 10^\{\-3\}\)Periodic Burgers 1D \(fixed viscosity\)\(8\.548×10−2\)±\(5\.901×10−3\)\\color\[rgb\]\{0,0,1\}\{\(8\.548\\times 10^\{\-2\}\)\\pm\(5\.901\\times 10^\{\-3\}\)\}\(9\.464×10−2\)±\(7\.872×10−3\)\(9\.464\\times 10^\{\-2\}\)\\pm\(7\.872\\times 10^\{\-3\}\)\(1\.356×10−1\)±\(8\.765×10−3\)\(1\.356\\times 10^\{\-1\}\)\\pm\(8\.765\\times 10^\{\-3\}\)Square\-pulse advection 1D\(2\.081×10−1\)±\(5\.286×10−3\)\\color\[rgb\]\{0,0,1\}\{\(2\.081\\times 10^\{\-1\}\)\\pm\(5\.286\\times 10^\{\-3\}\)\}\(2\.135×10−1\)±\(1\.751×10−3\)\(2\.135\\times 10^\{\-1\}\)\\pm\(1\.751\\times 10^\{\-3\}\)\(2\.406×10−1\)±\(4\.072×10−3\)\(2\.406\\times 10^\{\-1\}\)\\pm\(4\.072\\times 10^\{\-3\}\)Table 3:Five\-seed global relativeL2L\_\{2\}errors for nonlinear, shock, and Burgers benchmarks \(mean±\\pmsample standard deviation\)\. The best mean in each row is shown in blue\.BenchmarkFM\-OperatorFM\-DeepONet \(short\)DeepONet \(short\)FM\-DeepONetDeepONetShift\-DeepONetCubic conservation 1D0\.1656±0\.0133\\color\[rgb\]\{0,0,1\}\{0\.1656\\pm 0\.0133\}0\.1889±0\.02520\.1889\\pm 0\.02520\.2600±0\.01400\.2600\\pm 0\.01400\.1788±0\.00610\.1788\\pm 0\.00610\.2569±0\.01980\.2569\\pm 0\.01980\.2456±0\.02940\.2456\\pm 0\.0294Buckley–Leverett 1D0\.1245±0\.00730\.1245\\pm 0\.00730\.1100±0\.0062\\color\[rgb\]\{0,0,1\}\{0\.1100\\pm 0\.0062\}0\.1450±0\.00190\.1450\\pm 0\.00190\.1156±0\.00960\.1156\\pm 0\.00960\.1286±0\.00280\.1286\\pm 0\.00280\.1238±0\.00550\.1238\\pm 0\.0055Kuramoto–Sivashinsky 1D0\.0277±0\.00160\.0277\\pm 0\.00160\.0203±0\.00420\.0203\\pm 0\.00420\.0900±0\.01400\.0900\\pm 0\.01400\.0191±0\.0024\\color\[rgb\]\{0,0,1\}\{0\.0191\\pm 0\.0024\}0\.0207±0\.00460\.0207\\pm 0\.00460\.0835±0\.09120\.0835\\pm 0\.0912Square\-pulse advection 1D0\.2038±0\.0041\\color\[rgb\]\{0,0,1\}\{0\.2038\\pm 0\.0041\}0\.2102±0\.00730\.2102\\pm 0\.00730\.3249±0\.00640\.3249\\pm 0\.00640\.2216±0\.00890\.2216\\pm 0\.00890\.2414±0\.00400\.2414\\pm 0\.00400\.2274±0\.00270\.2274\\pm 0\.0027

Table 4:Global relativeL2L\_\{2\}errors over three independent seeds, reported as mean±\\pmsample standard deviation\. The best mean in each row is shown in blue\. All models use the same training and test functions within each benchmark\.
##### FM\-Operator vs\. Shift\-DeepONet\.

Shift\-DeepONet is designed for solutions with moving discontinuities or steep fronts\. Our Shift\-DeepONet implementation has a 64\-dimensional latent basis\. Its branch and trunk networks each use two hidden SiLU layers of width 128\. Two additional networks take the sensor values and generate an elementwise scale and shift in\(x,t\)\(x,t\)for every latent basis function, followed by a small coordinate\-dependent bias network\. This lets each basis function follow a front as its position changes\. The resulting model has about 141K parameters\. On the Burgers equation with Dirichlet boundary conditions andν=0\.001\\nu=0\.001, FM\-Operator uses about 43K parameters and obtains0\.0815±0\.00240\.0815\\pm 0\.0024, compared with0\.1132±0\.01570\.1132\\pm 0\.0157for Shift\-DeepONet\. This is a 28\.0% reduction in mean error, as shown in Table[5](https://arxiv.org/html/2607.28762#S5.T5)\.

ModelRelativeL2L\_\{2\}errorFM\-Operator0\.0815±0\.0024\\color\[rgb\]\{0,0,1\}\{0\.0815\\pm 0\.0024\}Shift\-DeepONet0\.1132±0\.01570\.1132\\pm 0\.0157Table 5:Comparison of FM\-Operator and Shift\-DeepONet on 1D Burgers with Dirichlet boundary conditions andν=0\.001\\nu=0\.001\. Results are global relativeL2L\_\{2\}errors over three independent seeds \(mean±\\pmsample standard deviation\)\.
##### Comparison on sharp\-gradient equations\.

Table[4](https://arxiv.org/html/2607.28762#S5.T4)gives a broader comparison on four equations with sharp gradients\. FM\-Operator remains competitive with only about 43K parameters: it gives the lowest error on cubic conservation and square\-pulse advection\. Shift\-DeepONet also performs well, particularly on Buckley–Leverett and square\-pulse advection\. The reduced DeepONet, with about 58K parameters, performs poorly on all four benchmarks\. Increasing DeepONet to about 108K parameters improves the results, but it is still worse than the best FM\-based model on cubic conservation, Buckley–Leverett, and square\-pulse advection\. On Kuramoto–Sivashinsky, the original DeepONet is competitive, while FM\-DeepONet gives the lowest error\. Figure[4](https://arxiv.org/html/2607.28762#S5.F4)provides an intuitive visualization of these results\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/four_benchmark_global_rel_l2.png)Figure 4:Global relativeL2L^\{2\}error comparison on four sharp\-gradient benchmark equations\.

## 6Discussion and Limitations

The results show two settings in which explicit feature interactions are particularly useful\. For PINNs, the advantage becomes clearer when the solution depends on several spatial coordinates\. FM\-PINN lowers the error on 16 of the 18 reported problems, with reductions of 43\.98% on 4D convection–diffusion, 55\.43% on 5D linear advection, and 40\.12% on 4D Burgers\. The weaker results on 2D Darcy and 2D reaction–diffusion also show that adding interaction terms does not improve every equation\.

For operator learning, the strongest improvements occur on conservation laws and other problems with shocks or steep fronts\. The sensor values describe the input function, while the query coordinates determine where that function is evaluated\. Modeling their interaction directly helps the network follow solution structures whose position and shape change with the input\. In contrast, splitting the query coordinates into separate fields changes the error only slightly on most benchmarks\. This suggests that sensor–query interactions matter more than interactions among the query coordinates themselves\.

The comparisons also show that the improvement is not simply caused by a larger network\. FM\-Operator uses about 43K parameters and outperforms the 108K\-parameter DeepONet on most nonlinear benchmarks in the five\-seed study\. The 58K\-parameter short DeepONet performs poorly on the four sharp\-gradient tests, whereas FM\-DeepONet improves the branch–trunk model by adding explicit interactions\. Shift\-DeepONet provides another strong baseline for moving fronts: it adapts its trunk basis through input\-dependent coordinate transformations, but on Burgers withν=0\.001\\nu=0\.001, FM\-Operator achieves a 28\.0% lower mean error with fewer parameters\.

The benefit is less consistent on smooth operator\-learning problems\. DeepONet is better on many Heat, Poisson, and Helmholtz settings, as shown in the supplementary results\. The current FM\-PINN comparisons are based on single runs, the main operator table uses five seeds, and the detailed model comparison uses three\. Future work should repeat the PINN experiments across seeds, test more PDE families, and study adaptive field grouping so that interaction terms are added only where they are useful\.

## 7Conclusion

This work introduces factorization\-machine\-style feature interactions to physics\-informed neural networks and neural operators\. FM\-PINN models interactions among spatial and temporal coordinates, while FM\-Operator and FM\-DeepONet extend the same idea to sensor values, query coordinates, summary statistics, and physical parameters\.

FM\-PINN improves accuracy on most of the reported multidimensional PDEs\. For neural operator learning, FM\-Operator is especially effective on nonlinear equations with shocks and steep gradients, despite using fewer parameters than DeepONet and Shift\-DeepONet\. FM\-DeepONet further shows that the interaction module can also strengthen the standard branch–trunk architecture\. The results on smooth equations are more mixed, indicating that explicit interactions are most valuable when the solution operator contains strong dependencies between the input function and the query location\. Overall, feature interaction modeling provides a simple and parameter\-efficient way to improve PDE learning in the problems where these dependencies are central\.

## References

- Y\. Chen, Y\. Lin, X\. Sun, C\. Yuan, and Z\. Gao \(2025\)Tensor decomposition\-based neural operator with dynamic mode decomposition for parameterized time\-dependent problems\.Journal of Computational Physics533,pp\. 113996\.External Links:ISSN 0021\-9991,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.jcp.2025.113996),[Link](https://www.sciencedirect.com/science/article/pii/S0021999125002797)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- J\. Cho, S\. Nam, H\. Yang, S\. Yun, Y\. Hong, and E\. Park \(2023\)Separable pinn: mitigating the curse of dimensionality in physics\-informed neural networks\.External Links:2211\.08761,[Link](https://arxiv.org/abs/2211.08761)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1),[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px1.p1.1)\.
- S\. Fazliani, Z\. Frangella, Y\. Chu, and M\. Udell \(2025\)Enhancing physics\-informed neural networks through feature engineering\.Transactions on Machine Learning Research\.Note:OpenReview: tS4vGObuCACited by:[§1](https://arxiv.org/html/2607.28762#S1.p3.1),[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- P\. Hadorn \(2022\)Shift\-deeponet: extending deep operator networks for discontinuous output functions\.External Links:[Link](https://api.semanticscholar.org/CorpusID:247862596)Cited by:[Appendix C](https://arxiv.org/html/2607.28762#A3.SSx2.SSSx7.p1.5)\.
- X\. He and T\. Chua \(2017\)Neural factorization machines for sparse predictive analytics\.InProceedings of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval,pp\. 355–364\.External Links:[Document](https://dx.doi.org/10.1145/3077136.3080777)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1),[§4](https://arxiv.org/html/2607.28762#S4.SSx1.p1.4)\.
- Z\. Hu, A\. D\. Jagtap, G\. E\. Karniadakis, and K\. Kawaguchi \(2023\)Augmented physics\-informed neural networks \(APINNs\): a gating network\-based soft domain decomposition methodology\.Engineering Applications of Artificial Intelligence126,pp\. 107183\.External Links:[Document](https://dx.doi.org/10.1016/j.engappai.2023.107183)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1),[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px1.p1.1)\.
- A\. D\. Jagtap and G\. E\. Karniadakis \(2020\)Extended physics\-informed neural networks \(XPINNs\): a generalized space\-time domain decomposition based deep learning framework for nonlinear partial differential equations\.Communications in Computational Physics28\(5\),pp\. 2002–2041\.External Links:[Document](https://dx.doi.org/10.4208/cicp.OA-2020-0164)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1),[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px1.p1.1)\.
- J\. Kossaifi, N\. Kovachki, K\. Azizzadenesheli, and A\. Anandkumar \(2023\)Multi\-grid tensorized fourier neural operator for high\-resolution pdes\.External Links:2310\.00120,[Link](https://arxiv.org/abs/2310.00120)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1),[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px2.p1.1)\.
- S\. Lanthaler, R\. Molinaro, P\. Hadorn, and S\. Mishra \(2023\)Nonlinear reconstruction for operator learning of PDEs with discontinuities\.InInternational Conference on Learning Representations,Cited by:[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px2.p1.1)\.
- Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar \(2021\)Fourier neural operator for parametric partial differential equations\.InInternational Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p2.1),[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- Z\. Li, H\. Zheng, N\. Kovachki, D\. Jin, H\. Chen, B\. Liu, K\. Azizzadenesheli, and A\. Anandkumar \(2023\)Physics\-informed neural operator for learning partial differential equations\.External Links:2111\.03794,[Link](https://arxiv.org/abs/2111.03794)Cited by:[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px2.p1.1)\.
- J\. Lian, X\. Zhou, F\. Zhang, Z\. Chen, X\. Xie, and G\. Sun \(2018\)xDeepFM: combining explicit and implicit feature interactions for recommender systems\.InProceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining,pp\. 1754–1763\.External Links:[Document](https://dx.doi.org/10.1145/3219819.3220023)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- Y\. Liu, Z\. Zhang, and H\. Schaeffer \(2023\)PROSE: predicting operators and symbolic expressions using multimodal transformers\.External Links:2309\.16816,[Link](https://arxiv.org/abs/2309.16816)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis \(2021\)Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.Nature Machine Intelligence3\(3\),pp\. 218–229\.External Links:[Document](https://dx.doi.org/10.1038/s42256-021-00302-5)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p2.1)\.
- B\. Moseley, A\. Markham, and T\. Nissen\-Meyer \(2023\)Finite basis physics\-informed neural networks \(fbpinns\): a scalable domain decomposition approach for solving differential equations\.Advances in Computational Mathematics49\(4\)\.External Links:ISSN 1572\-9044,[Link](http://dx.doi.org/10.1007/s10444-023-10065-9),[Document](https://dx.doi.org/10.1007/s10444-023-10065-9)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- M\. A\. Rahman, Z\. E\. Ross, and K\. Azizzadenesheli \(2023\)U\-no: u\-shaped neural operators\.External Links:2204\.11127,[Link](https://arxiv.org/abs/2204.11127)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis \(2019\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.10.045)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p1.1)\.
- S\. Rendle \(2010\)Factorization machines\.In2010 IEEE International Conference on Data Mining,pp\. 995–1000\.External Links:[Document](https://dx.doi.org/10.1109/ICDM.2010.127)Cited by:[Appendix B](https://arxiv.org/html/2607.28762#A2.p1.2),[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- M\. Tancik, P\. P\. Srinivasan, B\. Mildenhall, S\. Fridovich\-Keil, N\. Raghavan, U\. Singhal, R\. Ramamoorthi, J\. T\. Barron, and R\. Ng \(2020\)Fourier features let networks learn high frequency functions in low dimensional domains\.InAdvances in Neural Information Processing Systems,Vol\.33,pp\. 7537–7547\.Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p3.1),[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- A\. Tran, A\. Mathews, L\. Xie, and C\. S\. Ong \(2023\)Factorized fourier neural operators\.External Links:2111\.13802,[Link](https://arxiv.org/abs/2111.13802)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- S\. Wang, H\. Wang, and P\. Perdikaris \(2021\)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 Engineering384,pp\. 113938\.External Links:ISSN 0045\-7825,[Link](http://dx.doi.org/10.1016/j.cma.2021.113938),[Document](https://dx.doi.org/10.1016/j.cma.2021.113938)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p3.1),[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- J\. Xiao, H\. Ye, X\. He, H\. Zhang, F\. Wu, and T\. Chua \(2017\)Attentional factorization machines: learning the weight of feature interactions via attention networks\.InProceedings of the 26th International Joint Conference on Artificial Intelligence,pp\. 3119–3125\.External Links:[Document](https://dx.doi.org/10.24963/ijcai.2017/435)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1)\.
- M\. Zhu, J\. Sun, Z\. Zhang, H\. Schaeffer, and L\. Lu \(2025\)PI\-mfm: physics\-informed multimodal foundation model for solving partial differential equations\.External Links:2512\.23056,[Link](https://arxiv.org/abs/2512.23056)Cited by:[§1](https://arxiv.org/html/2607.28762#S1.p4.1),[§2](https://arxiv.org/html/2607.28762#S2.SS0.SSS0.Px2.p1.1)\.

## Appendix AMethodology

In this section, we introduce the three network architectures proposed in this paper: FM\-PINN, FM\-DeepONet, and FM\-Operator, which are illustrated in Figures[5](https://arxiv.org/html/2607.28762#A1.F5),[6](https://arxiv.org/html/2607.28762#A1.F6), and[7](https://arxiv.org/html/2607.28762#A1.F7), respectively\.

### FM\-PINN

![Refer to caption](https://arxiv.org/html/2607.28762v1/FMPINN.drawio.png)Figure 5:FM\-PINNFM\-PINN is a physics\-informed neural solver that factorizes the feature construction stage before applying the solution decoder\. Instead of concatenating all coordinate components into a single input vector, FM\-PINN treats each component as an individual field and embeds the fields separately before modeling their interactions\.

For a query pointyy, its input is represented as a collection ofMMfields:

ℱ​\(y\)=\{Fi​\(y\)\}i=1M,Fi​\(y\)∈ℝdi,\\mathcal\{F\}\(y\)=\\left\\\{F\_\{i\}\(y\)\\right\\\}\_\{i=1\}^\{M\},\\qquad F\_\{i\}\(y\)\\in\\mathbb\{R\}^\{d\_\{i\}\},whereFi​\(y\)F\_\{i\}\(y\)denotes the feature vector associated with theii\-th field anddid\_\{i\}is its input dimension\. For one\-dimensional time\-dependent equations, we set

F1​\(y\)=x,F2​\(y\)=tT,F\_\{1\}\(y\)=x,\\qquad F\_\{2\}\(y\)=\\frac\{t\}\{T\},whereTTdenotes the terminal time\. For two\-dimensional elliptic problems, we instead use

F1​\(y\)=x1,F2​\(y\)=x2\.F\_\{1\}\(y\)=x\_\{1\},\\qquad F\_\{2\}\(y\)=x\_\{2\}\.
Each field is independently mapped into a commondd\-dimensional embedding space through a field\-specific affine transformation:

ei​\(y\)=Wi​Fi​\(y\)\+βi,i=1,…,M,e\_\{i\}\(y\)=W\_\{i\}F\_\{i\}\(y\)\+\\beta\_\{i\},\\qquad i=1,\\ldots,M,\(1\)withWi∈ℝd×diW\_\{i\}\\in\\mathbb\{R\}^\{d\\times d\_\{i\}\},βi∈ℝd\\beta\_\{i\}\\in\\mathbb\{R\}^\{d\}, andei​\(y\)∈ℝde\_\{i\}\(y\)\\in\\mathbb\{R\}^\{d\}\. Here,ei​\(y\)e\_\{i\}\(y\)denotes the embedding vector of theii\-th field\. Mapping all fields to the same embedding dimension enables their subsequent feature\-interaction operations\.

The field embeddings are then passed to an NFM\-style bi\-interaction pooling layer,

zBI=12​\[\(∑i=1Mei\)⊙\(∑i=1Mei\)−∑i=1Mei⊙ei\]\.z\_\{\\mathrm\{BI\}\}=\\frac\{1\}\{2\}\\left\[\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\\odot\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\-\\sum\_\{i=1\}^\{M\}e\_\{i\}\\odot e\_\{i\}\\right\]\.\(2\)
The final prediction is obtained by applying an MLP decoder to the pooled interaction vector,

uθ​\(y\)=MLP​\(zBI\)\.u\_\{\\theta\}\(y\)=\\mathrm\{MLP\}\(z\_\{\\mathrm\{BI\}\}\)\.\(3\)Thus, FM\-PINN decouples continuous field embedding, explicit pairwise interactions, and nonlinear solution decoding\.

### FM\-DeepONet

![Refer to caption](https://arxiv.org/html/2607.28762v1/FM_DeepOnet.drawio.png)Figure 6:FM\-DeepONetFM\-DeepONet is an operator learning architecture that extends the standard DeepONet branch\-trunk framework by incorporating explicit factorization machine style feature interactions\. Concretely, it leverages an NFM\-based bi\-interaction module to capture interactions between branch and trunk features\.

The standard DeepONet prediction is given by the normalized inner product

𝒢θ​\(s\)​\(y\)=1p​⟨b,ttrunk⟩\+β\.\\mathcal\{G\}\_\{\\theta\}\(s\)\(y\)=\\frac\{1\}\{\\sqrt\{p\}\}\\langle b,t\_\{\\mathrm\{trunk\}\}\\rangle\+\\beta\.\(4\)
FM\-DeepONet preserves the branch\-trunk architecture of DeepONet while introducing an additional NFM interaction module that enhances the branch and trunk feature representations prior to the final inner\-product output\. Specifically, we construct a field set

ℱ=\{b,ttrunk,q​\(s\),η\},\\mathcal\{F\}=\\\{b,\\ t\_\{\\mathrm\{trunk\}\},\\ q\(s\),\\ \\eta\\\},whereq​\(s\)q\(s\)contains simple global statistics of the sensor vector,

q​\(s\)=\[mean​\(s\),std​\(s\),max⁡\|s\|\]∈ℝ3,q\(s\)=\\big\[\\mathrm\{mean\}\(s\),\\ \\mathrm\{std\}\(s\),\\ \\max\|s\|\\big\]\\in\\mathbb\{R\}^\{3\},andη\\etadenotes an optional scalar or low\-dimensional physical parameter, such as viscosity, wave speed, or diffusion coefficient\. If no such parameter is used, or if it is already included in the sensor input, the fieldη\\etais omitted\.

Each fieldFi∈ℱF\_\{i\}\\in\\mathcal\{F\}is mapped to a common interaction dimension via the same affine embedding as in Eq\. \([1](https://arxiv.org/html/2607.28762#A1.E1)\):

ei=Wi​Fi\+βi,ei∈ℝd\.e\_\{i\}=W\_\{i\}F\_\{i\}\+\\beta\_\{i\},\\qquad e\_\{i\}\\in\\mathbb\{R\}^\{d\}\.
Following the same NFM bi\-interaction pooling defined in Eq\. \([2](https://arxiv.org/html/2607.28762#A1.E2)\), the pooled interaction vector is

zBI=12​\[\(∑i=1Mei\)⊙\(∑i=1Mei\)−∑i=1Mei⊙ei\]\.z\_\{\\mathrm\{BI\}\}=\\frac\{1\}\{2\}\\left\[\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\\odot\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\-\\sum\_\{i=1\}^\{M\}e\_\{i\}\\odot e\_\{i\}\\right\]\.This vector explicitly aggregates pairwise interactions among the branch representation, trunk representation and sensor statistics\.

After that, the interaction vector is first processed by a multilayer perceptron \(MLP\) to obtain a refined representationrvectorr\_\{\\mathrm\{vector\}\}:

rvector=MLP​\(zBI\)\.r\_\{\\mathrm\{vector\}\}=\\mathrm\{MLP\}\(z\_\{\\mathrm\{BI\}\}\)\.Subsequently, this representation is projected into correction terms for the branch and trunk features:

Δ​b=Pb​\(rvector\),Δ​ttrunk=Ptrunk​\(rvector\),\\Delta b=P\_\{b\}\(r\_\{\\mathrm\{vector\}\}\),\\qquad\\Delta t\_\{\\mathrm\{trunk\}\}=P\_\{\\mathrm\{trunk\}\}\(r\_\{\\mathrm\{vector\}\}\),wherePbP\_\{b\}andPtrunkP\_\{\\mathrm\{trunk\}\}are trainable affine projections with outputs inℝp\\mathbb\{R\}^\{p\}\. The corrected branch and trunk features are

b~=b\+Δ​b,t~trunk=ttrunk\+Δ​ttrunk\.\\tilde\{b\}=b\+\\Delta b,\\qquad\\tilde\{t\}\_\{\\mathrm\{trunk\}\}=t\_\{\\mathrm\{trunk\}\}\+\\Delta t\_\{\\mathrm\{trunk\}\}\.The FM\-DeepONet prediction is computed using the corrected features while preserving the inner\-product structure of DeepONet:

𝒢θ′​\(s\)​\(y\)=1p​⟨b~,t~trunk⟩\+β\.\\mathcal\{G\}\_\{\\theta\}^\{\\mathrm\{\{\}^\{\\prime\}\}\}\(s\)\(y\)=\\frac\{1\}\{\\sqrt\{p\}\}\\langle\\tilde\{b\},\\tilde\{t\}\_\{\\mathrm\{trunk\}\}\\rangle\+\\beta\.\(5\)Therefore, FM\-DeepONet employs the NFM interaction module as a feature\-refinement mechanism before performing the branch–trunk inner product\. In this way, the model retains the decompositional structure of DeepONet while enabling pairwise interactions among input\-function features, query\-coordinate features, and optional equation parameters to jointly influence the final prediction\.

### FM\-Operator

![Refer to caption](https://arxiv.org/html/2607.28762v1/figures/FM_Operator.drawio.png)Figure 7:FM\-Operator architecture\.FM\-Operator does not preserve the branch\-trunk inner\-product structure of DeepONet\. Instead, it jointly feeds the input function, query coordinates, and relevant physical or parametric features into an FM/NFM module, which models their factorized feature interactions and directly predicts the solution value at the target point\.

Given the finite\-dimensional representationssof input function and a query coordinateyy, FM\-Operator constructs the field set

ℱ=\{s,y,q​\(s\),η\},\\mathcal\{F\}=\\\{s,\\ y,\\ q\(s\),\\ \\eta\\\},
Following the same convention as in FM\-DeepONet,η\\etais treated as an optional separate field and is omitted when the corresponding physical parameters are already included in the sensor inputss\. LetFi∈ℝdiF\_\{i\}\\in\\mathbb\{R\}^\{d\_\{i\}\}denote theii\-th field inℱ\\mathcal\{F\}\. Following the same embedding scheme as in Eq\. \([1](https://arxiv.org/html/2607.28762#A1.E1)\), each field is mapped viaei=Wi​Fi\+βie\_\{i\}=W\_\{i\}F\_\{i\}\+\\beta\_\{i\}withWi∈ℝd×diW\_\{i\}\\in\\mathbb\{R\}^\{d\\times d\_\{i\}\},βi∈ℝd\\beta\_\{i\}\\in\\mathbb\{R\}^\{d\}, andei∈ℝde\_\{i\}\\in\\mathbb\{R\}^\{d\}\. The NFM bi\-interaction pooling layer, defined in Eq\. \([2](https://arxiv.org/html/2607.28762#A1.E2)\), formulates second\-order field interactions as

zBI=12​\[\(∑i=1Mei\)⊙\(∑i=1Mei\)−∑i=1Mei⊙ei\],z\_\{\\mathrm\{BI\}\}=\\frac\{1\}\{2\}\\left\[\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\\odot\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\-\\sum\_\{i=1\}^\{M\}e\_\{i\}\\odot e\_\{i\}\\right\],where⊙\\odotdenotes Hadamard product\. The pooled interaction vector is then passed through a network to produce the interaction output\. In addition, a linear wide term is also applied to the concatenated original fields\.

Therefore, the final prediction is

𝒢θ′​\(s\)​\(y\)=MLP​\(zBI\)\+Linear​\(Concat​\(F1,…,FM\)\)\.\\mathcal\{G\}\_\{\\theta\}^\{\\mathrm\{\{\}^\{\\prime\}\}\}\(s\)\(y\)=\\mathrm\{MLP\}\(z\_\{\\mathrm\{BI\}\}\)\+\\mathrm\{Linear\}\(\\mathrm\{Concat\}\(F\_\{1\},\\ldots,F\_\{M\}\)\)\.\(6\)
Both FM\-DeepONet and FM\-Operator are trained in a supervised operator\-learning setting\. Given training PDE instances\{si\}i=1N\\\{s\_\{i\}\\\}\_\{i=1\}^\{N\}, query points\{yj\}j=1Q\\\{y\_\{j\}\\\}\_\{j=1\}^\{Q\}, and reference solution valuesun​\(yj\)u\_\{n\}\(y\_\{j\}\), the objective is

ℒθ=1N​Q​∑i=1N∑j=1Q\|𝒢θ′​\(si\)​\(yj\)−𝒢​\(si\)​\(yj\)\|2\.\\mathcal\{L\}\_\{\\theta\}=\\frac\{1\}\{NQ\}\\sum\_\{i=1\}^\{N\}\\sum\_\{j=1\}^\{Q\}\\left\|\\mathcal\{G\}\_\{\\theta\}^\{\\mathrm\{\{\}^\{\\prime\}\}\}\(s\_\{i\}\)\(y\_\{j\}\)\-\\mathcal\{G\}\(s\_\{i\}\)\(y\_\{j\}\)\\right\|^\{2\}\.\(7\)

## Appendix BTheoretical Analysis of Cross\-Field Interactions

A Factorization Machine \(FM\)\(Rendle[2010](https://arxiv.org/html/2607.28762#bib.bib9)\)models pairwise interactions as⟨vi,vj⟩​xi​xj\\langle v\_\{i\},v\_\{j\}\\rangle\\,x\_\{i\}x\_\{j\}, which is fundamentally a low\-rank factorization of the pairwise coefficient matrix\. From a functional decomposition perspective, this naturally corresponds to the pairwise components of a second\-order ANOVA expansion—the difference being that a full ANOVA additionally imposes orthogonality conditions and requires a reference distribution, neither of which we pursue here\. Since PDE solution maps are typically locally smooth, we adopt the more direct perspective of a local Taylor expansion: the cross\-field mixed derivativesδ​Fi⊤​Hi​j​δ​Fj\\delta F\_\{i\}^\{\\top\}H\_\{ij\}\\,\\delta F\_\{j\}that appear in the expansion are exactly the pairwise terms that FM’s bi\-interaction pooling is designed to explicitly construct\. This section develops this line of reasoning\. We first present the corresponding interpretation for FM\-PINN, then for FM\-Operator, and finally discuss the conditions under which the interpretation holds and the degenerate cases where it does not\.

### Local Cross\-Field Expansion

Split the input intoMMfields,

F=\(F1,…,FM\)∈𝒳:=ℝd1×⋯×ℝdM,F=\(F\_\{1\},\\ldots,F\_\{M\}\)\\in\\mathcal\{X\}:=\\mathbb\{R\}^\{d\_\{1\}\}\\times\\cdots\\times\\mathbb\{R\}^\{d\_\{M\}\},whereFi∈ℝdiF\_\{i\}\\in\\mathbb\{R\}^\{d\_\{i\}\}is fieldiianddid\_\{i\}is its dimension\. Letg:𝒳→ℝg:\\mathcal\{X\}\\to\\mathbb\{R\}be the scalar target mapping\. For FM\-PINN,ggis the solution value at a coordinate for one fixed PDE instance\. For FM\-Operator,ggis the finite\-dimensional target obtained after sensor discretization\.

Fix a reference inputF0=\(F10,…,FM0\)F^\{0\}=\(F\_\{1\}^\{0\},\\ldots,F\_\{M\}^\{0\}\)and define

δ​Fi:=Fi−Fi0,δ​F:=\(δ​F1,…,δ​FM\),‖δ​F‖2:=∑i=1M‖δ​Fi‖22\.\\delta F\_\{i\}:=F\_\{i\}\-F\_\{i\}^\{0\},\\qquad\\delta F:=\(\\delta F\_\{1\},\\ldots,\\delta F\_\{M\}\),\\qquad\\\|\\delta F\\\|^\{2\}:=\\sum\_\{i=1\}^\{M\}\\\|\\delta F\_\{i\}\\\|\_\{2\}^\{2\}\.Assumeg∈C3​\(𝒰\)g\\in C^\{3\}\(\\mathcal\{U\}\)on a convex neighborhood𝒰\\mathcal\{U\}containing the segment fromF0F^\{0\}toF0\+δ​FF^\{0\}\+\\delta F\. Then

g​\(F0\+δ​F\)=\\displaystyle g\(F^\{0\}\+\\delta F\)=\{\}g​\(F0\)\+∑i=1M∇Fig​\(F0\)⊤​δ​Fi\\displaystyle g\(F^\{0\}\)\+\\sum\_\{i=1\}^\{M\}\\nabla\_\{F\_\{i\}\}g\(F^\{0\}\)^\{\\top\}\\delta F\_\{i\}\+12​∑i=1Mδ​Fi⊤​Hi​i​\(F0\)​δ​Fi\+∑1≤i<j≤Mδ​Fi⊤​Hi​j​\(F0\)​δ​Fj\+R3​\(δ​F\),\\displaystyle\+\\frac\{1\}\{2\}\\sum\_\{i=1\}^\{M\}\\delta F\_\{i\}^\{\\top\}H\_\{ii\}\(F^\{0\}\)\\delta F\_\{i\}\+\\sum\_\{1\\leq i<j\\leq M\}\\delta F\_\{i\}^\{\\top\}H\_\{ij\}\(F^\{0\}\)\\delta F\_\{j\}\+R\_\{3\}\(\\delta F\),\(8\)where

Hi​j​\(F0\):=∇Fi​Fj2g​\(F0\)∈ℝdi×djH\_\{ij\}\(F^\{0\}\):=\\nabla^\{2\}\_\{F\_\{i\}F\_\{j\}\}g\(F^\{0\}\)\\in\\mathbb\{R\}^\{d\_\{i\}\\times d\_\{j\}\}is the\(i,j\)\(i,j\)Hessian block\. IfsupF∈𝒰‖D3​g​\(F\)‖op≤C3\\sup\_\{F\\in\\mathcal\{U\}\}\\\|D^\{3\}g\(F\)\\\|\_\{\\mathrm\{op\}\}\\leq C\_\{3\}, withD3​gD^\{3\}gthe third derivative and∥⋅∥op\\\|\\cdot\\\|\_\{\\mathrm\{op\}\}its operator norm, then

\|R3​\(δ​F\)\|≤C36​‖δ​F‖3\.\|R\_\{3\}\(\\delta F\)\|\\leq\\frac\{C\_\{3\}\}\{6\}\\\|\\delta F\\\|^\{3\}\.
The cross\-field part is

Qcross​\(δ​F\):=∑1≤i<j≤Mδ​Fi⊤​Hi​j​\(F0\)​δ​Fj\.Q\_\{\\mathrm\{cross\}\}\(\\delta F\):=\\sum\_\{1\\leq i<j\\leq M\}\\delta F\_\{i\}^\{\\top\}H\_\{ij\}\(F^\{0\}\)\\delta F\_\{j\}\.\(9\)

### FM\-PINN

For a fixed PDE instance, FM\-PINN learns the coordinate\-to\-solution mapy↦u​\(y\)y\\mapsto u\(y\), whereyyis a spatial or space\-time coordinate\. For a steadyDD\-dimensional problem, treat each coordinate as one scalar field:Fi=xi∈ℝF\_\{i\}=x\_\{i\}\\in\\mathbb\{R\}fori=1,…,Di=1,\\ldots,D\. Withδ​xi:=xi−xi0\\delta x\_\{i\}:=x\_\{i\}\-x\_\{i\}^\{0\}, equation \([8](https://arxiv.org/html/2607.28762#A2.E8)\) becomes

u​\(x0\+δ​x\)=\\displaystyle u\(x^\{0\}\+\\delta x\)=\{\}u​\(x0\)\+∑i=1Duxi​\(x0\)​δ​xi\+12​∑i=1Duxi​xi​\(x0\)​\(δ​xi\)2\\displaystyle u\(x^\{0\}\)\+\\sum\_\{i=1\}^\{D\}u\_\{x\_\{i\}\}\(x^\{0\}\)\\delta x\_\{i\}\+\\frac\{1\}\{2\}\\sum\_\{i=1\}^\{D\}u\_\{x\_\{i\}x\_\{i\}\}\(x^\{0\}\)\(\\delta x\_\{i\}\)^\{2\}\+∑i<juxi​xj​\(x0\)​δ​xi​δ​xj\+R3\.\\displaystyle\+\\sum\_\{i<j\}u\_\{x\_\{i\}x\_\{j\}\}\(x^\{0\}\)\\delta x\_\{i\}\\delta x\_\{j\}\+R\_\{3\}\.\(10\)The cross\-coordinate structure is∑i<juxi​xj​\(x0\)​δ​xi​δ​xj\\sum\_\{i<j\}u\_\{x\_\{i\}x\_\{j\}\}\(x^\{0\}\)\\delta x\_\{i\}\\delta x\_\{j\}\. ADD\-dimensional input yields\(D2\)\\binom\{D\}\{2\}pairwise terms; a benefit is plausible only when the corresponding mixed derivatives are non\-negligible\.

For a one\-dimensional time\-dependent problem, use the two fieldsF1=xF\_\{1\}=xandF2=tF\_\{2\}=t\. The mixed term reduces toux​t​δ​x​δ​tu\_\{xt\}\\,\\delta x\\,\\delta t\.

### FM\-Operator

After sensor discretization, write the target as

g​\(s,y,η\)=𝒢​\(s,η\)​\(y\),g\(s,y,\\eta\)=\\mathcal\{G\}\(s,\\eta\)\(y\),wheres∈ℝnss\\in\\mathbb\{R\}^\{n\_\{s\}\}is the sensor vector,y∈ℝnyy\\in\\mathbb\{R\}^\{n\_\{y\}\}is the query coordinate,η∈ℝnη\\eta\\in\\mathbb\{R\}^\{n\_\{\\eta\}\}is an optional physical parameter, and𝒢\\mathcal\{G\}is the resulting finite\-dimensional solution operator\. Around\(s0,y0,η0\)\(s^\{0\},y^\{0\},\\eta^\{0\}\), define

δ​s:=s−s0,δ​y:=y−y0,δ​η:=η−η0\.\\delta s:=s\-s^\{0\},\\qquad\\delta y:=y\-y^\{0\},\\qquad\\delta\\eta:=\\eta\-\\eta^\{0\}\.Under the same localC3C^\{3\}condition, letg0g\_\{0\}be the value at the reference point,gs,gy,gηg\_\{s\},g\_\{y\},g\_\{\\eta\}the three gradients, andHp​qH\_\{pq\}the Hessian block withp,q∈\{s,y,η\}p,q\\in\\\{s,y,\\eta\\\}\. The local expansion is

g​\(s0\+δ​s,y0\+δ​y,η0\+δ​η\)=\\displaystyle g\(s^\{0\}\+\\delta s,y^\{0\}\+\\delta y,\\eta^\{0\}\+\\delta\\eta\)=\{\}g0\+gs⊤​δ​s\+gy⊤​δ​y\+gη⊤​δ​η\\displaystyle g\_\{0\}\+g\_\{s\}^\{\\top\}\\delta s\+g\_\{y\}^\{\\top\}\\delta y\+g\_\{\\eta\}^\{\\top\}\\delta\\eta\+12​δ​s⊤​Hs​s​δ​s\+12​δ​y⊤​Hy​y​δ​y\+12​δ​η⊤​Hη​η​δ​η\\displaystyle\+\\frac\{1\}\{2\}\\delta s^\{\\top\}H\_\{ss\}\\delta s\+\\frac\{1\}\{2\}\\delta y^\{\\top\}H\_\{yy\}\\delta y\+\\frac\{1\}\{2\}\\delta\\eta^\{\\top\}H\_\{\\eta\\eta\}\\delta\\eta\+δ​s⊤​Hs​y​δ​y\+δ​s⊤​Hs​η​δ​η\+δ​y⊤​Hy​η​δ​η\+R3,\\displaystyle\+\\delta s^\{\\top\}H\_\{sy\}\\delta y\+\\delta s^\{\\top\}H\_\{s\\eta\}\\delta\\eta\+\\delta y^\{\\top\}H\_\{y\\eta\}\\delta\\eta\+R\_\{3\},\(11\)where\|R3\|≤C36​\(‖δ​s‖22\+‖δ​y‖22\+‖δ​η‖22\)3/2\|R\_\{3\}\|\\leq\\frac\{C\_\{3\}\}\{6\}\(\\\|\\delta s\\\|\_\{2\}^\{2\}\+\\\|\\delta y\\\|\_\{2\}^\{2\}\+\\\|\\delta\\eta\\\|\_\{2\}^\{2\}\)^\{3/2\}under the usual third\-derivative bound\.

FM\-Operator treatsssas one field andyyas another\. Its main explicit cross\-field term is therefore

δ​s⊤​Hs​y​δ​y\.\\delta s^\{\\top\}H\_\{sy\}\\delta y\.Ifη\\etais a separate field, parameter\-related interactions are also available; ifη\\etais concatenated withss, it belongs to the sensor block instead\.

Linearity of the PDE operator inssdoes not remove the sensor–query interaction\. For example, if

𝒢​\(s\)​\(y\)=∑m=1nssm​ϕm​\(y\),\\mathcal\{G\}\(s\)\(y\)=\\sum\_\{m=1\}^\{n\_\{s\}\}s\_\{m\}\\phi\_\{m\}\(y\),withϕm\\phi\_\{m\}differentiable, then

∇s​s2g​\(s,y\)=0,∂2g∂sm​∂yk=∂ϕm​\(y\)∂yk≠0\.\\nabla^\{2\}\_\{ss\}g\(s,y\)=0,\\qquad\\frac\{\\partial^\{2\}g\}\{\\partial s\_\{m\}\\partial y\_\{k\}\}=\\frac\{\\partial\\phi\_\{m\}\(y\)\}\{\\partial y\_\{k\}\}\\neq 0\.Thus the relevant question is not whether the PDE is linear, but whether the joint map\(s,y\)↦𝒢​\(s\)​\(y\)\(s,y\)\\mapsto\\mathcal\{G\}\(s\)\(y\)depends substantially on both fields\.

### What Bi\-Interaction Pooling Represents

Map fieldiito an embedding

ei=Wi​Fi\+βi∈ℝr,e\_\{i\}=W\_\{i\}F\_\{i\}\+\\beta\_\{i\}\\in\\mathbb\{R\}^\{r\},whereWi∈ℝr×diW\_\{i\}\\in\\mathbb\{R\}^\{r\\times d\_\{i\}\},βi∈ℝr\\beta\_\{i\}\\in\\mathbb\{R\}^\{r\}, andrris the embedding dimension\. NFM pooling computes

zBI:=12​\[\(∑i=1Mei\)⊙\(∑i=1Mei\)−∑i=1Mei⊙ei\],z\_\{\\mathrm\{BI\}\}:=\\frac\{1\}\{2\}\\left\[\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\\odot\\left\(\\sum\_\{i=1\}^\{M\}e\_\{i\}\\right\)\-\\sum\_\{i=1\}^\{M\}e\_\{i\}\\odot e\_\{i\}\\right\],where⊙\\odotdenotes element\-wise multiplication\. Expanding the square yields

zBI=∑1≤i<j≤Mei⊙ej,z\_\{\\mathrm\{BI\}\}=\\sum\_\{1\\leq i<j\\leq M\}e\_\{i\}\\odot e\_\{j\},so the layer explicitly forms a product for every pair of distinct fields\.

For a scalar linear readoutf​\(F\)=a⊤​zBI\+bf\(F\)=a^\{\\top\}z\_\{\\mathrm\{BI\}\}\+bwitha∈ℝra\\in\\mathbb\{R\}^\{r\},b∈ℝb\\in\\mathbb\{R\}, letDa:=diag⁡\(a\)D\_\{a\}:=\\operatorname\{diag\}\(a\)\. The purely bilinear part for pairi<ji<jis

Fi⊤​Bi​j​Fj,Bi​j:=Wi⊤​Da​Wj,rank⁡\(Bi​j\)≤r,F\_\{i\}^\{\\top\}B\_\{ij\}F\_\{j\},\\qquad B\_\{ij\}:=W\_\{i\}^\{\\top\}D\_\{a\}W\_\{j\},\\qquad\\operatorname\{rank\}\(B\_\{ij\}\)\\leq r,sinceBi​jB\_\{ij\}factors throughℝr\\mathbb\{R\}^\{r\}\. Biases contribute linear and constant terms\. With a nonlinear decoder, the interaction vector is no longer a single bilinear form, but the pairwise structure remains explicit\.

experimental results serve as the primary evidence\.

## Appendix CExperiments

### FM\-PINN

In this section, we compare the performance of FM\-PINN and the standard PINN on several smooth high\-dimensional problems\. The goal is to demonstrate that FM\-PINN can improve the performance of PINNs on high\-dimensional problems by explicitly constructing second\-order feature interactions\.

#### Experiment Configuration

We compare a standard physics\-informed neural network, denoted as PINN, with the proposed feature\-interaction model, denoted as FM\-PINN\. In the main experiments, PINN uses a coordinate\-based MLP feature projection, while FM\-PINN uses a neural factorization\-machine \(NFM\) feature interaction layer\. All comparisons are conducted under the same balanced loss setting, using identical sampling strategies, optimization settings, training steps, and evaluation metrics\.

Table 6:Architectural difference between PINN and FM\-PINN\.ModelFeature representationPINNRaw coordinate projectionFM\-PINNNFM feature interaction
Note\.The shared decoder consists of Linear64→12864\\rightarrow 128, SiLU activation, four hidden layers with width 128 and SiLU activations, and a final Linear128→1128\\rightarrow 1output layer\.

The model hyperparameters are

feature dimension=64,hidden dimension=128,decoder depth=4\.\\text\{feature dimension\}=64,\\qquad\\text\{hidden dimension\}=128,\\qquad\\text\{decoder depth\}=4\.
All models are trained for30,00030\{,\}000optimization steps using AdamW with

γlr=10−3,weight decay=10−6\.\\gamma\_\{\\mathrm\{lr\}\}=10^\{\-3\},\\qquad\\text\{weight decay\}=10^\{\-6\}\.
Collocation points and boundary or initial points are sampled using Latin Hypercube Sampling\. The sampling configuration is

Nf=10000,Nbc/ic=400\.N\_\{f\}=10000,\\qquad N\_\{\\mathrm\{bc/ic\}\}=400\.
The mini\-batch sizes are

Bf=1024,Bbc=400,Bic=400\.B\_\{f\}=1024,\\qquad B\_\{\\mathrm\{bc\}\}=400,\\qquad B\_\{\\mathrm\{ic\}\}=400\.

#### Experiment Results

In the implementation of FM\-PINN, we further investigated the effect of the first\-order main\-effect term in NFM on model performance\. A standard NFM typically consists of a first\-order linear term,w0\+∑iwi​Fiw\_\{0\}\+\\sum\_\{i\}w\_\{i\}F\_\{i\}, together with higher\-order feature interaction terms\. We compared three settings: PINN, NFM\-PINN, and NFM\-PINN with the first\-order term\. As shown in Table[7](https://arxiv.org/html/2607.28762#A3.T7), retaining the first\-order term does not lead to consistent performance improvement, and its effect depends on the specific PDE type and dimensionality\. In particular, for several high\-dimensional problems, NFM\-PINN without the first\-order linear term achieves lower relativeL2L^\{2\}errors, suggesting a more stable empirical advantage\. Therefore, the subsequent analysis is mainly based on the interaction\-only NFM\-PINN, namely the variant without the first\-order linear term\.

Table 7:Ablation study of the first\-order term in NFM\-PINN\.PDEPINNNFM\-PINNNFM\-PINN with first\-order termHeat \(1D\)3\.135×10−043\.135\\text\{\\times\}\{10\}^\{\-04\}2\.806×10−042\.806\\text\{\\times\}\{10\}^\{\-04\}6\.560×10−046\.560\\text\{\\times\}\{10\}^\{\-04\}Poisson \(2D\)1\.850×10−021\.850\\text\{\\times\}\{10\}^\{\-02\}1\.782×10−021\.782\\text\{\\times\}\{10\}^\{\-02\}1\.413×10−021\.413\\text\{\\times\}\{10\}^\{\-02\}Wave \(1D\)1\.436×10−021\.436\\text\{\\times\}\{10\}^\{\-02\}1\.437×10−021\.437\\text\{\\times\}\{10\}^\{\-02\}1\.580×10−021\.580\\text\{\\times\}\{10\}^\{\-02\}Darcy \(2D\)2\.411×10−032\.411\\text\{\\times\}\{10\}^\{\-03\}3\.156×10−033\.156\\text\{\\times\}\{10\}^\{\-03\}2\.950×10−032\.950\\text\{\\times\}\{10\}^\{\-03\}Burgers \(1D,ν=0\.01\\nu=0\.01\)1\.249×10−021\.249\\text\{\\times\}\{10\}^\{\-02\}1\.237×10−021\.237\\text\{\\times\}\{10\}^\{\-02\}1\.233×10−021\.233\\text\{\\times\}\{10\}^\{\-02\}Heat \(4D\)8\.018×10−038\.018\\text\{\\times\}\{10\}^\{\-03\}5\.014×10−035\.014\\text\{\\times\}\{10\}^\{\-03\}7\.238×10−037\.238\\text\{\\times\}\{10\}^\{\-03\}Poisson \(4D\)1\.352×10−011\.352\\text\{\\times\}\{10\}^\{\-01\}1\.298×10−011\.298\\text\{\\times\}\{10\}^\{\-01\}1\.355×10−011\.355\\text\{\\times\}\{10\}^\{\-01\}Poisson \(10D\)4\.589×10−014\.589\\text\{\\times\}\{10\}^\{\-01\}4\.225×10−014\.225\\text\{\\times\}\{10\}^\{\-01\}4\.122×10−014\.122\\text\{\\times\}\{10\}^\{\-01\}Convection–Diffusion \(1D\)6\.718×10−046\.718\\text\{\\times\}\{10\}^\{\-04\}9\.985×10−049\.985\\text\{\\times\}\{10\}^\{\-04\}9\.862×10−049\.862\\text\{\\times\}\{10\}^\{\-04\}Convection–Diffusion \(2D\)1\.898×10−031\.898\\text\{\\times\}\{10\}^\{\-03\}1\.031×10−031\.031\\text\{\\times\}\{10\}^\{\-03\}1\.204×10−031\.204\\text\{\\times\}\{10\}^\{\-03\}Convection–Diffusion \(4D\)2\.602×10−022\.602\\text\{\\times\}\{10\}^\{\-02\}1\.458×10−021\.458\\text\{\\times\}\{10\}^\{\-02\}1\.488×10−021\.488\\text\{\\times\}\{10\}^\{\-02\}Convection–Diffusion \(5D\)3\.817×10−023\.817\\text\{\\times\}\{10\}^\{\-02\}3\.012×10−023\.012\\text\{\\times\}\{10\}^\{\-02\}3\.399×10−023\.399\\text\{\\times\}\{10\}^\{\-02\}Reaction–Diffusion \(1D\)6\.854×10−056\.854\\text\{\\times\}\{10\}^\{\-05\}2\.710×10−042\.710\\text\{\\times\}\{10\}^\{\-04\}2\.371×10−042\.371\\text\{\\times\}\{10\}^\{\-04\}Reaction–Diffusion \(2D\)4\.608×10−044\.608\\text\{\\times\}\{10\}^\{\-04\}5\.672×10−045\.672\\text\{\\times\}\{10\}^\{\-04\}3\.242×10−043\.242\\text\{\\times\}\{10\}^\{\-04\}Reaction–Diffusion \(4D\)1\.498×10−021\.498\\text\{\\times\}\{10\}^\{\-02\}1\.272×10−021\.272\\text\{\\times\}\{10\}^\{\-02\}1\.266×10−021\.266\\text\{\\times\}\{10\}^\{\-02\}Reaction–Diffusion \(5D\)3\.270×10−023\.270\\text\{\\times\}\{10\}^\{\-02\}2\.632×10−022\.632\\text\{\\times\}\{10\}^\{\-02\}3\.128×10−023\.128\\text\{\\times\}\{10\}^\{\-02\}Linear Advection \(1D\)7\.456×10−057\.456\\text\{\\times\}\{10\}^\{\-05\}6\.836×10−056\.836\\text\{\\times\}\{10\}^\{\-05\}1\.239×10−041\.239\\text\{\\times\}\{10\}^\{\-04\}Linear Advection \(2D\)4\.213×10−034\.213\\text\{\\times\}\{10\}^\{\-03\}2\.495×10−032\.495\\text\{\\times\}\{10\}^\{\-03\}6\.213×10−046\.213\\text\{\\times\}\{10\}^\{\-04\}Linear Advection \(4D\)2\.111×10−022\.111\\text\{\\times\}\{10\}^\{\-02\}7\.755×10−037\.755\\text\{\\times\}\{10\}^\{\-03\}8\.740×10−038\.740\\text\{\\times\}\{10\}^\{\-03\}Linear Advection \(5D\)3\.133×10−023\.133\\text\{\\times\}\{10\}^\{\-02\}1\.396×10−021\.396\\text\{\\times\}\{10\}^\{\-02\}1\.948×10−021\.948\\text\{\\times\}\{10\}^\{\-02\}Darcy \(4D\)1\.066×10−011\.066\\text\{\\times\}\{10\}^\{\-01\}1\.005×10−011\.005\\text\{\\times\}\{10\}^\{\-01\}1\.197×10−011\.197\\text\{\\times\}\{10\}^\{\-01\}Wave \(4D\)4\.988×10−014\.988\\text\{\\times\}\{10\}^\{\-01\}3\.149×10−013\.149\\text\{\\times\}\{10\}^\{\-01\}3\.032×10−013\.032\\text\{\\times\}\{10\}^\{\-01\}Burgers \(4D, pure coordinates\)6\.489×10−036\.489\\text\{\\times\}\{10\}^\{\-03\}3\.886×10−033\.886\\text\{\\times\}\{10\}^\{\-03\}3\.798×10−033\.798\\text\{\\times\}\{10\}^\{\-03\}##### Smooth Steady\-State Equations

We consider convection\-diffusion, reaction\-diffusion, and linear advection equations on

Ω=\[0,1\]d,d∈\{1,2,4,5\}\.\\Omega=\[0,1\]^\{d\},\\qquad d\\in\\\{1,2,4,5\\\}\.
The boundary condition is

u​\(x\)=u∗​\(x\),x∈∂Ω\.u\(x\)=u^\{\\ast\}\(x\),\\qquad x\\in\\partial\\Omega\.
The manufactured exact solution is

u∗​\(x\)=∑m=15cm​∏j=1dsin⁡\(km,j​π​xj\)\.u^\{\\ast\}\(x\)=\\sum\_\{m=1\}^\{5\}c\_\{m\}\\prod\_\{j=1\}^\{d\}\\sin\(k\_\{m,j\}\\pi x\_\{j\}\)\.
The coefficients are

\(c1,…,c5\)=\(1\.00,0\.35,−0\.25,0\.20,0\.10\)\.\(c\_\{1\},\\ldots,c\_\{5\}\)=\(1\.00,\\ 0\.35,\\ \-0\.25,\\ 0\.20,\\ 0\.10\)\.
The frequency vectors are truncated to dimensiondd:

k1=\(1,1,…,1\),k2=\(2,1,…,1\),k3=\(1,2,1,…,1\),k\_\{1\}=\(1,1,\\ldots,1\),\\quad k\_\{2\}=\(2,1,\\ldots,1\),\\quad k\_\{3\}=\(1,2,1,\\ldots,1\),
k4=\(2,1,2,1,…,1\),k5=\(3,1,2,1,…,1\)\.k\_\{4\}=\(2,1,2,1,\\ldots,1\),\\quad k\_\{5\}=\(3,1,2,1,\\ldots,1\)\.- •For convection\-diffusion Equation, we solve −ν​Δ​u\+β⋅∇u=f​\(x\),x∈Ω,\-\\nu\\Delta u\+\\beta\\cdot\\nabla u=f\(x\),\\qquad x\\in\\Omega, with The velocity vector is truncated from β=\(1\.0,−0\.7,0\.5,−0\.3,0\.2\)\.\\beta=\(1\.0,\\ \-0\.7,\\ 0\.5,\\ \-0\.3,\\ 0\.2\)\. The forcing term is f​\(x\)=−ν​Δ​u∗​\(x\)\+β⋅∇u∗​\(x\)\.f\(x\)=\-\\nu\\Delta u^\{\\ast\}\(x\)\+\\beta\\cdot\\nabla u^\{\\ast\}\(x\)\.
- •For reaction\-diffusion Equation, we solve −ν​Δ​u\+ρ​\(u3−u\)=f​\(x\),x∈Ω,\-\\nu\\Delta u\+\\rho\(u^\{3\}\-u\)=f\(x\),\\qquad x\\in\\Omega, with ν=0\.05,ρ=1\.0\.\\nu=0\.05,\\qquad\\rho=1\.0\. The forcing term is f​\(x\)=−ν​Δ​u∗​\(x\)\+ρ​\(\(u∗​\(x\)\)3−u∗​\(x\)\)\.f\(x\)=\-\\nu\\Delta u^\{\\ast\}\(x\)\+\\rho\\left\(\(u^\{\\ast\}\(x\)\)^\{3\}\-u^\{\\ast\}\(x\)\\right\)\.
- •For linear advection Equation, we solve β⋅∇u=f​\(x\),x∈Ω,\\beta\\cdot\\nabla u=f\(x\),\\qquad x\\in\\Omega, with f​\(x\)=β⋅∇u∗​\(x\)\.f\(x\)=\\beta\\cdot\\nabla u^\{\\ast\}\(x\)\.

##### Observations\.

NFM\-PINN shows no advantage over PINN in one dimension: for 1D convection–diffusion theL2L^\{2\}error is 48\.62% higher, and for 1D linear advection it improves by only 8\.31%\. A single coordinate field provides no pairs for the bi\-interaction layer, which degenerates to a constant embedding \(Section 4, single\-field degeneracy\)\.

The benefit emerges in higher dimensions\. For 2D convection–diffusion and linear advection, NFM\-PINN reduces theL2L^\{2\}error by 45\.69% and 40\.78%, respectively\. For 4D and 5D linear advection, the reductions reach 63\.26% and 55\.43%\. The number of available cross\-coordinate pairs grows from 1 \(2D\) to 6 \(4D\) and 10 \(5D\)\.

Heatmap and 3D\-slice visualizations in Figures[10](https://arxiv.org/html/2607.28762#A3.F10)–[13](https://arxiv.org/html/2607.28762#A3.F13)provide qualitative comparisons\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_convection_diffusion1d.png)\(a\)Convection Diffusion 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_convection_diffusion2d.png)\(b\)Convection Diffusion 2D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_convection_diffusion4d.png)\(c\)Convection Diffusion 4D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_convection_diffusion5d.png)\(d\)Convection\-diffusion 5D

Figure 8:Comparison between PINN and FM\-PINN on Convection Diffusion\.![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_linear_advection1d.png)\(a\)Linear Advection 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_linear_advection2d.png)\(b\)Linear Advection 2D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_linear_advection4d.png)\(c\)Linear Advection 4D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/04_balanced_smooth_linear_advection5d.png)\(d\)Linear Advection 5D

Figure 9:Comparison between PINN and FM\-PINN on Linear Advection\.![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/other/convection_diffusion2d_mlp_pinn_best_triptych.png)\(a\)PINN Convection Diffusion 2D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/other/convection_diffusion2d_nfm_pinn_best_triptych.png)\(b\)FM\-PINN Convection Diffusion 2D

Figure 10:Comparison between PINN and FM\-PINN on Convection Diffusion 2D\.![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/other/linear_advection2d_mlp_pinn_best_triptych.png)\(a\)PINN Linear Advection
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/other/linear_advection2d_nfm_pinn_best_triptych.png)\(b\)FM\-PINN Linear Advection

Figure 11:Comparison between PINN and FM\-PINN on Linear Advection\.![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/4d/convection_diffusion4d_mlp_pinn_best_volume.png)\(a\)PINN Convection Diffusion 4D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/4d/convection_diffusion4d_nfm_pinn_best_volume.png)\(b\)FM\-PINN Convection Diffusion 4D

Figure 12:Comparison between PINN and FM\-PINN on Convection Diffusion 4D\.![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/4d/linear_advection4d_mlp_pinn_best_volume.png)\(a\)PINN Linear Advection 4D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_smooth/heatmap/4d/linear_advection4d_nfm_pinn_best_volume.png)\(b\)FM\-PINN CLinear Advection 4D

Figure 13:Comparison between PINN and FM\-PINN on Linear Advection 4D\.
##### Heat Equations

For the one\-dimensional heat equation,

x∈\[0,1\],t∈\[0,0\.25\],x\\in\[0,1\],\\qquad t\\in\[0,0\.25\],
we solve

ut−ν​ux​x=0,ν=0\.05\.u\_\{t\}\-\\nu u\_\{xx\}=0,\\qquad\\nu=0\.05\.
The exact solution is

u∗​\(x,t\)=∑m=15cm​sin⁡\(km​π​x\)​exp⁡\(−ν​\(km​π\)2​t\),u^\{\\ast\}\(x,t\)=\\sum\_\{m=1\}^\{5\}c\_\{m\}\\sin\(k\_\{m\}\\pi x\)\\exp\\left\(\-\\nu\(k\_\{m\}\\pi\)^\{2\}t\\right\),
where

\(km,cm\)=\(1,1\.00\),\(2,−0\.35\),\(3,0\.25\),\(4,−0\.15\),\(5,0\.08\)\.\(k\_\{m\},c\_\{m\}\)=\(1,1\.00\),\\ \(2,\-0\.35\),\\ \(3,0\.25\),\\ \(4,\-0\.15\),\\ \(5,0\.08\)\.
The initial condition is

u​\(x,0\)=u∗​\(x,0\),u\(x,0\)=u^\{\\ast\}\(x,0\),
and the boundary condition is

u​\(0,t\)=u​\(1,t\)=0\.u\(0,t\)=u\(1,t\)=0\.
For the four\-dimensional heat equation,

\(x,y,z\)∈\[0,1\]3,t∈\[0,0\.25\],\(x,y,z\)\\in\[0,1\]^\{3\},\\qquad t\\in\[0,0\.25\],
we solve

ut−ν​Δ​u=0\.u\_\{t\}\-\\nu\\Delta u=0\.
The exact solution is

u∗​\(x,y,z,t\)=∑mcm​sin⁡\(km,1​π​x\)​sin⁡\(km,2​π​y\)​sin⁡\(km,3​π​z\)​exp⁡\(−ν​λm​t\),u^\{\\ast\}\(x,y,z,t\)=\\sum\_\{m\}c\_\{m\}\\sin\(k\_\{m,1\}\\pi x\)\\sin\(k\_\{m,2\}\\pi y\)\\sin\(k\_\{m,3\}\\pi z\)\\exp\(\-\\nu\\lambda\_\{m\}t\),
where

λm=π2​\(km,12\+km,22\+km,32\)\.\\lambda\_\{m\}=\\pi^\{2\}\(k\_\{m,1\}^\{2\}\+k\_\{m,2\}^\{2\}\+k\_\{m,3\}^\{2\}\)\.
Initial and boundary conditions are imposed fromu∗u^\{\\ast\}\.

##### Observations\.

FM\-PINN achieves anL2L^\{2\}error of5\.014×10−35\.014\\times 10^\{\-3\}on the 4D heat equation, a 37\.46% reduction from PINN’s8\.018×10−38\.018\\times 10^\{\-3\}\(Figure[14](https://arxiv.org/html/2607.28762#A3.F14)\)\. The four coordinate fields\(x,y,z,t/T\)\(x,y,z,t/T\)yield\(42\)=6\\binom\{4\}\{2\}=6pairwise interactions\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_hbw/02_balanced_main_heat1d.png)\(a\)Heat 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_hbw/02_balanced_main_heat4d.png)\(b\)Heat 4D

Figure 14:Comparison between PINN and FM\-PINN on Heat\.
##### Burgers Equations

For the one\-dimensional Burgers equation,

x∈\[−1,1\],t∈\[0,1\],x\\in\[\-1,1\],\\qquad t\\in\[0,1\],
we solve

ut\+u​ux−ν​ux​x=0,ν=0\.01\.u\_\{t\}\+uu\_\{x\}\-\\nu u\_\{xx\}=0,\\qquad\\nu=0\.01\.
The initial condition is

u​\(x,0\)=−sin⁡\(π​x\),u\(x,0\)=\-\\sin\(\\pi x\),
and the boundary condition is

u​\(−1,t\)=u​\(1,t\)=0\.u\(\-1,t\)=u\(1,t\)=0\.
For the four\-dimensional Burgers\-type equation,

\(x,y,z\)∈\[0,1\]3,t∈\[0,1\],\(x,y,z\)\\in\[0,1\]^\{3\},\\qquad t\\in\[0,1\],
we solve

ut\+u​\(ux\+uy\+uz\)−ν​Δ​u=f​\(x,y,z,t\),ν=0\.01\.u\_\{t\}\+u\(u\_\{x\}\+u\_\{y\}\+u\_\{z\}\)\-\\nu\\Delta u=f\(x,y,z,t\),\\qquad\\nu=0\.01\.
The manufactured solution is

u∗​\(x,y,z,t\)=∑mcm​sin⁡\(km,1​π​x\)​sin⁡\(km,2​π​y\)​sin⁡\(km,3​π​z\)​exp⁡\(−ν​λm​t\),u^\{\\ast\}\(x,y,z,t\)=\\sum\_\{m\}c\_\{m\}\\sin\(k\_\{m,1\}\\pi x\)\\sin\(k\_\{m,2\}\\pi y\)\\sin\(k\_\{m,3\}\\pi z\)\\exp\(\-\\nu\\lambda\_\{m\}t\),
where

λm=π2​\(km,12\+km,22\+km,32\)\.\\lambda\_\{m\}=\\pi^\{2\}\(k\_\{m,1\}^\{2\}\+k\_\{m,2\}^\{2\}\+k\_\{m,3\}^\{2\}\)\.
The terms are

\(km,cm\)=\(\(1,1,1\),−0\.90\),\(\(2,1,1\),0\.25\),\(\(1,2,1\),−0\.20\),\(\(1,1,2\),0\.15\)\.\(k\_\{m\},c\_\{m\}\)=\(\(1,1,1\),\-0\.90\),\(\(2,1,1\),0\.25\),\(\(1,2,1\),\-0\.20\),\(\(1,1,2\),0\.15\)\.
The forcing term is

f=ut∗\+u∗​\(ux∗\+uy∗\+uz∗\)−ν​Δ​u∗\.f=u\_\{t\}^\{\\ast\}\+u^\{\\ast\}\(u\_\{x\}^\{\\ast\}\+u\_\{y\}^\{\\ast\}\+u\_\{z\}^\{\\ast\}\)\-\\nu\\Delta u^\{\\ast\}\.
Initial and boundary conditions are imposed fromu∗u^\{\\ast\}\.

##### Observations\.

FM\-PINN and PINN perform nearly identically on the 1D viscous Burgers problem with a shock\-like steep gradient \(Figure[15](https://arxiv.org/html/2607.28762#A3.F15)\)\. Cross\-coordinate interactions do not resolve the inherent difficulty of residual\-based training in steep\-gradient regimes: the residual becomes highly localized, and the training error is dominated by collocation coverage, derivative accuracy, and spectral bias\.

For the 4D forced Burgers equation with a smooth manufactured solution, FM\-PINN reduces theL2L^\{2\}error by 40\.12%\. The contrast highlights that the interaction mechanism helps when the solution depends smoothly on multiple coordinates but does not compensate for the optimization challenges of sharp layers\. In a separate data\-assisted experiment \(not shown\), adding a supervised data termℒdata\\mathcal\{L\}\_\{\\mathrm\{data\}\}restores the advantage of FM\-PINN on the 1D shock problem, suggesting that the limitation lies in the residual\-based objective rather than the architecture itself\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_hbw/02_balanced_main_burgers1d_nu0.01.png)\(a\)Burgers 1D \- shock
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_hbw/03_balanced_extra4d_burgers4d_pure.png)\(b\)Burgers 4D \-smooth

Figure 15:Comparison between PINN and FM\-PINN on Burgers\.
##### Wave Equations

For the one\-dimensional wave equation,

x∈\[0,1\],t∈\[0,1\],x\\in\[0,1\],\\qquad t\\in\[0,1\],
we solve

ut​t−cwave2​ux​x=0,cwave=1\.u\_\{tt\}\-c\_\{\\mathrm\{wave\}\}^\{2\}u\_\{xx\}=0,\\qquad c\_\{\\mathrm\{wave\}\}=1\.
The exact solution is

u∗​\(x,t\)=sin⁡\(π​x\)​cos⁡\(π​t\)\+0\.5​sin⁡\(2​π​x\)​cos⁡\(2​π​t\)\.u^\{\\ast\}\(x,t\)=\\sin\(\\pi x\)\\cos\(\\pi t\)\+0\.5\\sin\(2\\pi x\)\\cos\(2\\pi t\)\.
The initial displacement is

u​\(x,0\)=u∗​\(x,0\),u\(x,0\)=u^\{\\ast\}\(x,0\),
the initial velocity is

and the boundary condition is

u​\(0,t\)=u​\(1,t\)=0\.u\(0,t\)=u\(1,t\)=0\.
For the four\-dimensional wave equation,

\(x,y,z\)∈\[0,1\]3,t∈\[0,1\],\(x,y,z\)\\in\[0,1\]^\{3\},\\qquad t\\in\[0,1\],
we solve

ut​t−cwave2​Δ​u=0,cwave=1\.u\_\{tt\}\-c\_\{\\mathrm\{wave\}\}^\{2\}\\Delta u=0,\\qquad c\_\{\\mathrm\{wave\}\}=1\.
The exact solution is

u∗​\(x,y,z,t\)=∑mcm​sin⁡\(km,1​π​x\)​sin⁡\(km,2​π​y\)​sin⁡\(km,3​π​z\)​cos⁡\(ωm​t\),u^\{\\ast\}\(x,y,z,t\)=\\sum\_\{m\}c\_\{m\}\\sin\(k\_\{m,1\}\\pi x\)\\sin\(k\_\{m,2\}\\pi y\)\\sin\(k\_\{m,3\}\\pi z\)\\cos\(\\omega\_\{m\}t\),
where

ωm=c​π​km,12\+km,22\+km,32\.\\omega\_\{m\}=c\\pi\\sqrt\{k\_\{m,1\}^\{2\}\+k\_\{m,2\}^\{2\}\+k\_\{m,3\}^\{2\}\}\.
The terms are

\(km,cm\)=\(\(1,1,1\),1\.00\),\(\(2,1,1\),0\.35\),\(\(1,2,3\),−0\.25\),\(\(2,2,1\),0\.20\),\(\(3,1,2\),0\.10\)\.\(k\_\{m\},c\_\{m\}\)=\(\(1,1,1\),1\.00\),\(\(2,1,1\),0\.35\),\(\(1,2,3\),\-0\.25\),\(\(2,2,1\),0\.20\),\(\(3,1,2\),0\.10\)\.
The initial displacement, initial velocity, and boundary conditions are imposed fromu∗u^\{\\ast\}\. In particular,

ut​\(x,y,z,0\)=0\.u\_\{t\}\(x,y,z,0\)=0\.

##### Observations\.

FM\-PINN achieves a 36\.87%L2L^\{2\}error reduction on the 4D wave equation \(3\.149×10−13\.149\\times 10^\{\-1\}vs\. PINN’s4\.988×10−14\.988\\times 10^\{\-1\}; Figure[16](https://arxiv.org/html/2607.28762#A3.F16)\), while performing slightly worse in 1D\. The 4D problem has four coordinate fields and six pairwise interactions; the 1D time\-dependent problem has only two fields with one interaction pair\. The large absolute errors of both models indicate that the 4D wave equation remains challenging\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_hbw/02_balanced_main_wave1d.png)\(a\)Wave 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/picture_fmpinn_hbw/03_balanced_extra4d_wave4d.png)\(b\)Wave 4D

Figure 16:Comparison between PINN and FM\-PINN on Wave\.

#### Result Analysis

Table[8](https://arxiv.org/html/2607.28762#A3.T8)summarizes the error reductions across all benchmark problems\. Three patterns emerge\.

FM\-PINN improves high\-dimensional smooth problems\.On 4D and 5D convection–diffusion, linear advection, and reaction–diffusion, FM\-PINN reduces theL2L^\{2\}error by 30–63% relative to PINN\. The 4D forced Burgers and 4D heat equations show similar gains \(40\.12% and 37\.46%\)\. As the number of coordinate fields grows from 2 to 5, the available cross\-coordinate pairs increase from 1 to 10\.

Limited benefit in 1D\.On 1D steady\-state problems, FM\-PINN provides no advantage \(single\-field degeneracy\)\. On 1D time\-dependent problems, a single\(x,t\)\(x,t\)interaction pair yields inconsistent improvement\.

No improvement on shock\-like solutions\.For the 1D Burgers equation with a steep gradient, FM\-PINN and PINN perform identically under the residual\-based objective\. The interaction mechanism does not address the optimization difficulty posed by highly localized residuals\. A data\-assisted variant \(not shown\) recovers the advantage of FM\-PINN on this problem\.

Overall, FM\-PINN is best understood as an architectural enhancement for high\-dimensional smooth PDEs where cross\-coordinate dependence is substantial, rather than a remedy for the optimization challenges of shock\-like solutions\.

Moreover, in the data\-assisted setting, FM\-PINN becomes more effective after introducing the supervised termℒdata\\mathcal\{L\}\_\{\\mathrm\{data\}\}\. This suggests that direct observations in the steep\-gradient region can alleviate the residual\-optimization bottleneck and enable the interaction architecture to use its representational capacity more effectively\.

We also incorporated AFM\- and xDeepFM\-style interaction modules into the PINN framework\. In the experiments considered here, AFM exhibits unstable and strongly problem\-dependent optimization behavior\. Although it obtains competitive accuracy on a small number of equations, it generally underperforms the NFM\-based variant\. The xDeepFM\-based PINN introduces a deeper explicit interaction structure, but this additional complexity substantially increases the difficulty of optimizing the physics\-informed residual and frequently results in slow or unsuccessful convergence\. Among the investigated interaction modules, NFM therefore provides the most favorable empirical balance between structured second\-order interaction modeling and optimization stability\.

Table 8:RelativeL2L^\{2\}error comparison between PINN and FM\-PINN\. The error reduction is computed as\(PINN−FM​\-​PINN\)/PINN×100%\(\\mathrm\{PINN\}\-\\mathrm\{FM\}\\mbox\{\-\}\\mathrm\{PINN\}\)/\\mathrm\{PINN\}\\times 100\\%; negative values indicate that FM\-PINN has a higher error than PINN\.PDEPINNFM\-PINNError reductionPINN/FM\-PINNResultHeat \(1D\)3\.135×10−043\.135\\text\{\\times\}\{10\}^\{\-04\}2\.806×10−042\.806\\text\{\\times\}\{10\}^\{\-04\}10\.49%1\.121\.12FM\-PINN lower errorPoisson \(2D\)1\.850×10−021\.850\\text\{\\times\}\{10\}^\{\-02\}1\.782×10−021\.782\\text\{\\times\}\{10\}^\{\-02\}3\.64%1\.041\.04FM\-PINN lower errorWave \(1D\)1\.436×10−021\.436\\text\{\\times\}\{10\}^\{\-02\}1\.437×10−021\.437\\text\{\\times\}\{10\}^\{\-02\}\-0\.06%1\.001\.00FM\-PINN higher errorDarcy \(2D\)2\.411×10−032\.411\\text\{\\times\}\{10\}^\{\-03\}3\.156×10−033\.156\\text\{\\times\}\{10\}^\{\-03\}\-30\.91%0\.760\.76FM\-PINN higher errorBurgers \(1D,ν=0\.01\\nu=0\.01\)1\.249×10−021\.249\\text\{\\times\}\{10\}^\{\-02\}1\.237×10−021\.237\\text\{\\times\}\{10\}^\{\-02\}0\.96%1\.011\.01FM\-PINN lower errorHeat \(4D\)8\.018×10−038\.018\\text\{\\times\}\{10\}^\{\-03\}5\.014×10−035\.014\\text\{\\times\}\{10\}^\{\-03\}37\.46%1\.601\.60FM\-PINN lower errorPoisson \(4D\)1\.352×10−011\.352\\text\{\\times\}\{10\}^\{\-01\}1\.298×10−011\.298\\text\{\\times\}\{10\}^\{\-01\}3\.99%1\.041\.04FM\-PINN lower errorPoisson \(10D\)4\.589×10−014\.589\\text\{\\times\}\{10\}^\{\-01\}4\.225×10−014\.225\\text\{\\times\}\{10\}^\{\-01\}7\.92%1\.091\.09FM\-PINN lower errorConvection–Diffusion \(1D\)6\.718×10−046\.718\\text\{\\times\}\{10\}^\{\-04\}9\.985×10−049\.985\\text\{\\times\}\{10\}^\{\-04\}\-48\.62%0\.670\.67FM\-PINN higher errorConvection–Diffusion \(2D\)1\.898×10−031\.898\\text\{\\times\}\{10\}^\{\-03\}1\.031×10−031\.031\\text\{\\times\}\{10\}^\{\-03\}45\.69%1\.841\.84FM\-PINN lower errorConvection–Diffusion \(4D\)2\.602×10−022\.602\\text\{\\times\}\{10\}^\{\-02\}1\.458×10−021\.458\\text\{\\times\}\{10\}^\{\-02\}43\.98%1\.791\.79FM\-PINN lower errorConvection–Diffusion \(5D\)3\.817×10−023\.817\\text\{\\times\}\{10\}^\{\-02\}3\.012×10−023\.012\\text\{\\times\}\{10\}^\{\-02\}21\.08%1\.271\.27FM\-PINN lower errorReaction–Diffusion \(1D\)6\.854×10−056\.854\\text\{\\times\}\{10\}^\{\-05\}2\.710×10−042\.710\\text\{\\times\}\{10\}^\{\-04\}\-295\.37%0\.250\.25FM\-PINN higher errorReaction–Diffusion \(2D\)4\.608×10−044\.608\\text\{\\times\}\{10\}^\{\-04\}5\.672×10−045\.672\\text\{\\times\}\{10\}^\{\-04\}\-23\.09%0\.810\.81FM\-PINN higher errorReaction–Diffusion \(4D\)1\.498×10−021\.498\\text\{\\times\}\{10\}^\{\-02\}1\.272×10−021\.272\\text\{\\times\}\{10\}^\{\-02\}15\.14%1\.181\.18FM\-PINN lower errorReaction–Diffusion \(5D\)3\.270×10−023\.270\\text\{\\times\}\{10\}^\{\-02\}2\.632×10−022\.632\\text\{\\times\}\{10\}^\{\-02\}19\.51%1\.241\.24FM\-PINN lower errorLinear Advection \(1D\)7\.456×10−057\.456\\text\{\\times\}\{10\}^\{\-05\}6\.836×10−056\.836\\text\{\\times\}\{10\}^\{\-05\}8\.31%1\.091\.09FM\-PINN lower errorLinear Advection \(2D\)4\.213×10−034\.213\\text\{\\times\}\{10\}^\{\-03\}2\.495×10−032\.495\\text\{\\times\}\{10\}^\{\-03\}40\.78%1\.691\.69FM\-PINN lower errorLinear Advection \(4D\)2\.111×10−022\.111\\text\{\\times\}\{10\}^\{\-02\}7\.755×10−037\.755\\text\{\\times\}\{10\}^\{\-03\}63\.26%2\.722\.72FM\-PINN lower errorLinear Advection \(5D\)3\.133×10−023\.133\\text\{\\times\}\{10\}^\{\-02\}1\.396×10−021\.396\\text\{\\times\}\{10\}^\{\-02\}55\.43%2\.242\.24FM\-PINN lower errorDarcy \(4D\)1\.066×10−011\.066\\text\{\\times\}\{10\}^\{\-01\}1\.005×10−011\.005\\text\{\\times\}\{10\}^\{\-01\}5\.73%1\.061\.06FM\-PINN lower errorWave \(4D\)4\.988×10−014\.988\\text\{\\times\}\{10\}^\{\-01\}3\.149×10−013\.149\\text\{\\times\}\{10\}^\{\-01\}36\.87%1\.581\.58FM\-PINN lower errorBurgers \(4D\)6\.489×10−036\.489\\text\{\\times\}\{10\}^\{\-03\}3\.886×10−033\.886\\text\{\\times\}\{10\}^\{\-03\}40\.12%1\.671\.67FM\-PINN lower error

### FM\-Operator and FM\-DeepONet

We evaluate FM\-DeepONet and FM\-Operator against DeepONet on a suite of nonlinear conservation laws and smooth PDEs\. All models use identical training data, optimization, and evaluation protocols\.

#### Experiment Configuration

Unless otherwise stated, each operator\-learning dataset contains

Ntrain=256,Nval=64,Ntest=64N\_\{\\mathrm\{train\}\}=256,\\qquad N\_\{\\mathrm\{val\}\}=64,\\qquad N\_\{\\mathrm\{test\}\}=64PDE instances or input functions\. For non\-parametric operator benchmarks, each input function is represented byns=64n\_\{s\}=64sensor values\. For the parametric Burgers experiment, the operator input additionally contains the normalized viscosity and two Dirichlet boundary values\.

For DeepONet and FM\-DeepONet, the branch input dimension is6464in non\-parametric benchmarks and6767in the parametric Burgers experiment\.

Each instance is evaluated on a full space\-time grid with

nx=128,nt=80,n\_\{x\}=128,\\qquad n\_\{t\}=80,which provides128×80128\\times 80supervised query values per instance\.

For DeepONet and FM\-DeepONet, the branch and trunk networks use hidden width128128, latent dimensionp=128p=128, and depth33\. For FM\-DeepONet, the NFM interaction embedding dimension is128128, and the interaction hidden width is128128\. For FM\-Operator, each raw field is embedded into a common128128\-dimensional interaction space, followed by an NFM interaction decoder with hidden width128128\. For xDeepFM\-based ablation variants, the CIN module uses88compressed interaction fields\. The activation function is SiLU\.

All models are trained for30,00030\{,\}000optimization steps using AdamW with

γlr=10−3,weight decay=10−6\.\\gamma\_\{\\mathrm\{lr\}\}=10^\{\-3\},\\qquad\\text\{weight decay\}=10^\{\-6\}\.Each mini\-batch samples

Bfunc=32B\_\{\\mathrm\{func\}\}=32training functions and

Bpoint=128B\_\{\\mathrm\{point\}\}=128query points per function, resulting in40964096supervised pointwise samples per optimization step\.

Table 9:Number of trainable parameters for the operator\-learning models\.EquationDeepONetFM\-DeepONet\(XDeepFM\)FM\-DeepONet\(NFM\)FM\-Operator\(XDeepFM\)FM\-Operator\(NFM\)Cubic conservation107,777453,401207,361339,09742,439Buckley–Leverett107,777453,401207,361339,09742,439Square\-pulse advection107,777453,401207,361339,09742,439Fixed\-viscosity periodic Burgers107,777453,777207,617339,72942,696

#### Experiment Results

Table 10:Best\-train\-loss test MSE and global relativeL2L\_\{2\}for the corrected operator experiments\. Each cell reports test MSE on the first line and global relativeL2L\_\{2\}on the second line\.EquationDeepONetFM\-DeepONet\(XDeepFM\)FM\-DeepONet\(NFM\)FM\-Operator\(XDeepFM\)FM\-Operator\(NFM\)Periodic BurgersMSE:1\.643×10−31\.643\\times 10^\{\-3\}L2L\_\{2\}:1\.282×10−11\.282\\times 10^\{\-1\}MSE:8\.986×10−48\.986\\times 10^\{\-4\}L2L\_\{2\}:9\.479×10−29\.479\\times 10^\{\-2\}MSE:2\.672×𝟏𝟎−𝟒\\mathbf\{2\.672\\times 10^\{\-4\}\}L2L\_\{2\}:5\.169×𝟏𝟎−𝟐\\mathbf\{5\.169\\times 10^\{\-2\}\}MSE:4\.029×10−34\.029\\times 10^\{\-3\}L2L\_\{2\}:2\.007×10−12\.007\\times 10^\{\-1\}MSE:7\.116×10−47\.116\\times 10^\{\-4\}L2L\_\{2\}:8\.436×10−28\.436\\times 10^\{\-2\}Cubic conservationMSE:7\.713×10−37\.713\\times 10^\{\-3\}L2L\_\{2\}:2\.319×10−12\.319\\times 10^\{\-1\}MSE:7\.060×10−37\.060\\times 10^\{\-3\}L2L\_\{2\}:2\.218×10−12\.218\\times 10^\{\-1\}MSE:2\.865×𝟏𝟎−𝟑\\mathbf\{2\.865\\times 10^\{\-3\}\}L2L\_\{2\}:1\.413×𝟏𝟎−𝟏\\mathbf\{1\.413\\times 10^\{\-1\}\}MSE:1\.356×10−21\.356\\times 10^\{\-2\}L2L\_\{2\}:3\.074×10−13\.074\\times 10^\{\-1\}MSE:3\.459×10−33\.459\\times 10^\{\-3\}L2L\_\{2\}:1\.553×10−11\.553\\times 10^\{\-1\}Buckley–LeverettMSE:2\.049×10−32\.049\\times 10^\{\-3\}L2L\_\{2\}:8\.649×10−28\.649\\times 10^\{\-2\}MSE:2\.131×10−32\.131\\times 10^\{\-3\}L2L\_\{2\}:8\.819×10−28\.819\\times 10^\{\-2\}MSE:1\.018×𝟏𝟎−𝟑\\mathbf\{1\.018\\times 10^\{\-3\}\}L2L\_\{2\}:6\.097×𝟏𝟎−𝟐\\mathbf\{6\.097\\times 10^\{\-2\}\}MSE:2\.160×10−32\.160\\times 10^\{\-3\}L2L\_\{2\}:8\.881×10−28\.881\\times 10^\{\-2\}MSE:1\.285×10−31\.285\\times 10^\{\-3\}L2L\_\{2\}:6\.848×10−26\.848\\times 10^\{\-2\}Square\-pulse advectionMSE:1\.270×10−21\.270\\times 10^\{\-2\}L2L\_\{2\}:2\.404×10−12\.404\\times 10^\{\-1\}MSE:1\.046×10−21\.046\\times 10^\{\-2\}L2L\_\{2\}:2\.182×10−12\.182\\times 10^\{\-1\}MSE:1\.075×10−21\.075\\times 10^\{\-2\}L2L\_\{2\}:2\.212×10−12\.212\\times 10^\{\-1\}MSE:1\.071×10−21\.071\\times 10^\{\-2\}L2L\_\{2\}:2\.207×10−12\.207\\times 10^\{\-1\}MSE:9\.846×𝟏𝟎−𝟑\\mathbf\{9\.846\\times 10^\{\-3\}\}L2L\_\{2\}:2\.117×𝟏𝟎−𝟏\\mathbf\{2\.117\\times 10^\{\-1\}\}Kuramoto–SivashinskyMSE:2\.362×𝟏𝟎−𝟓\\mathbf\{2\.362\\times 10^\{\-5\}\}L2L\_\{2\}:1\.583×𝟏𝟎−𝟐\\mathbf\{1\.583\\times 10^\{\-2\}\}MSE:3\.506×10−43\.506\\times 10^\{\-4\}L2L\_\{2\}:6\.100×10−26\.100\\times 10^\{\-2\}MSE:4\.016×10−54\.016\\times 10^\{\-5\}L2L\_\{2\}:2\.065×10−22\.065\\times 10^\{\-2\}MSE:2\.391×10−32\.391\\times 10^\{\-3\}L2L\_\{2\}:1\.593×10−11\.593\\times 10^\{\-1\}MSE:9\.160×10−59\.160\\times 10^\{\-5\}L2L\_\{2\}:3\.118×10−23\.118\\times 10^\{\-2\}LWR traffic flowMSE:2\.152×10−32\.152\\times 10^\{\-3\}L2L\_\{2\}:8\.251×10−28\.251\\times 10^\{\-2\}MSE:9\.943×𝟏𝟎−𝟒\\mathbf\{9\.943\\times 10^\{\-4\}\}L2L\_\{2\}:5\.609×𝟏𝟎−𝟐\\mathbf\{5\.609\\times 10^\{\-2\}\}MSE:1\.027×10−31\.027\\times 10^\{\-3\}L2L\_\{2\}:5\.700×10−25\.700\\times 10^\{\-2\}MSE:1\.142×10−31\.142\\times 10^\{\-3\}L2L\_\{2\}:6\.010×10−26\.010\\times 10^\{\-2\}MSE:1\.048×10−31\.048\\times 10^\{\-3\}L2L\_\{2\}:5\.758×10−25\.758\\times 10^\{\-2\}

![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/total_bar.png)Figure 17:relative global relativeL2L\_\{2\}error##### Cubic Conservation Law

The cubic conservation experiment considers a nonlinear scalar conservation law with small diffusion,

ut\+\(u33\)x=ν​ux​x,x∈\[0,1\),t∈\[0,1\],u\_\{t\}\+\\left\(\\frac\{u^\{3\}\}\{3\}\\right\)\_\{x\}=\\nu u\_\{xx\},\\qquad x\\in\[0,1\),\\quad t\\in\[0,1\],with

The data are generated on a periodic grid\. The numerical fluxes and the diffusion stencil use wrap\-around indexing, corresponding to the periodic identification

u​\(0,t\)=u​\(1,t\)\.u\(0,t\)=u\(1,t\)\.
The initial functions are smooth periodic random fields,

u​\(x,0\)=α\+∑k=16\(ak​sin⁡\(2​π​k​x\)\+bk​cos⁡\(2​π​k​x\)\),u\(x,0\)=\\alpha\\;\+\\;\\sum\_\{k=1\}^\{6\}\\left\(a\_\{k\}\\sin\(2\\pi kx\)\+b\_\{k\}\\cos\(2\\pi kx\)\\right\),where the coefficients decay with1/k1/k\. In this experimentα=0\\alpha=0, and the initial fields are normalized and scaled with amplitude0\.90\.9\. Reference solutions are generated by a periodic finite\-volume method using a Rusanov numerical flux, with an additional periodic diffusion term\.

Results\.See Table[10](https://arxiv.org/html/2607.28762#A3.T10)and Figure[17](https://arxiv.org/html/2607.28762#A3.F17)for quantitative comparisons, and Figure[18](https://arxiv.org/html/2607.28762#A3.F18)for solution profiles\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/cubic_conservation1d_deeponet_best_trainloss_test0.png)\(a\)DeepONet Cubic Conservation Law
![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/cubic_conservation1d_nfm_deeponet_best_trainloss_test0.png)\(b\)FM\-Operator Cubic Conservation Law
![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/cubic_conservation1d_nfm_feature_deeponet_best_trainloss_test0.png)\(c\)FM\-DeepONet\(NFM\) Cubic Conservation Law

Figure 18:Comparison between DeepONet, FM\-Operator and FM\-DeepONet\(NFM\) on Cubic Conservation Law\.
##### Buckley–Leverett Equation

The Buckley–Leverett experiment solves

ut\+f​\(u\)x=0,x∈\[0,1\),t∈\[0,0\.6\],u\_\{t\}\+f\(u\)\_\{x\}=0,\\qquad x\\in\[0,1\),\\quad t\\in\[0,0\.6\],on a periodic grid\. The flux is

f​\(u\)=u2u2\+a​\(1−u\)2,a=0\.5\.f\(u\)=\\frac\{u^\{2\}\}\{u^\{2\}\+a\(1\-u\)^\{2\}\},\\qquad a=0\.5\.The periodic boundary is implemented through wrap\-around finite\-volume fluxes, equivalently identifying

u​\(0,t\)=u​\(1,t\)\.u\(0,t\)=u\(1,t\)\.
The initial condition is a smooth periodic random field with offset0\.50\.5and amplitude0\.420\.42, clipped to the interval\[0\.02,0\.98\]\[0\.02,0\.98\]\. The reference solution is generated by the same periodic finite\-volume procedure with Rusanov flux and no artificial diffusion\. This problem contains nonlinear wave steepening and is used to test the ability of the operator models to represent shock\-like solution features\.

Results\.See Table[10](https://arxiv.org/html/2607.28762#A3.T10)and Figure[19](https://arxiv.org/html/2607.28762#A3.F19)\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/buckley_leverett1d_deeponet_best_trainloss_test0.png)\(a\)DeepONet Buckley–Leverett Equation
![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/buckley_leverett1d_nfm_deeponet_best_trainloss_test0.png)\(b\)FM\-Operator Buckley–Leverett Equation
![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/buckley_leverett1d_nfm_feature_deeponet_best_trainloss_test0.png)\(c\)FM\-DeepONet\(NFM\) Buckley–Leverett Equation

Figure 19:Comparison between DeepONet, FM\-Operator and FM\-DeepONet\(NFM\) on Buckley–Leverett Equation\.
##### Square\-Pulse Linear Advection

The square\-pulse advection experiment considers the one\-dimensional linear advection equation

ut\+c​ux=0,x∈\[0,1\),t∈\[0,1\],u\_\{t\}\+cu\_\{x\}=0,\\qquad x\\in\[0,1\),\\quad t\\in\[0,1\],on a periodic domain, with constant advection speed

The periodic boundary condition is imposed as

u​\(0,t\)=u​\(1,t\)\.u\(0,t\)=u\(1,t\)\.
For each PDE instance, the input function is an initial square pulse

u​\(x,0\)=A​1\{\|\(\(x−xc\+1/2\)mod1\)−1/2\|≤w2\},u\(x,0\)=A\\,\\mathbf\{1\}\_\{\\left\\\{\\left\|\\left\(\(x\-x\_\{c\}\+1/2\)\\bmod 1\\right\)\-1/2\\right\|\\leq\\frac\{w\}\{2\}\\right\\\}\},where the amplitude, width, and initial center are sampled independently as

A∼𝒰​\(0\.6,1\.4\),w∼𝒰​\(0\.12,0\.28\),xc∼𝒰​\(0\.20,0\.80\)\.A\\sim\\mathcal\{U\}\(0\.6,1\.4\),\\qquad w\\sim\\mathcal\{U\}\(0\.12,0\.28\),\\qquad x\_\{c\}\\sim\\mathcal\{U\}\(0\.20,0\.80\)\.The reference solution is generated analytically by translating the initial pulse with periodic wrap\-around\. Specifically, the pulse center evolves as

xc​\(t\)=\(xc\+c​t\)mod1,x\_\{c\}\(t\)=\(x\_\{c\}\+ct\)\\bmod 1,and the solution is

u​\(x,t\)=A​1\{\|\(\(x−xc​\(t\)\+1/2\)mod1\)−1/2\|≤w2\}\.u\(x,t\)=A\\,\\mathbf\{1\}\_\{\\left\\\{\\left\|\\left\(\(x\-x\_\{c\}\(t\)\+1/2\)\\bmod 1\\right\)\-1/2\\right\|\\leq\\frac\{w\}\{2\}\\right\\\}\}\.This setting provides a nonsmooth transport benchmark for evaluating whether the learned operator can propagate discontinuous input functions under periodic advection\.

Results\.See Table[10](https://arxiv.org/html/2607.28762#A3.T10)and Figure[20](https://arxiv.org/html/2607.28762#A3.F20)\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/square_advection1d_deeponet_best_trainloss_test0.png)\(a\)DeepONet Square linear advection 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/square_advection1d_nfm_deeponet_best_trainloss_test0.png)\(b\)FM\-Operator Square linear advection 1D

Figure 20:Comparison between DeepONet and FM\-Operator on Square linear advection 1D \.
##### Fixed\-Viscosity Periodic Burgers Equation

The fixed\-viscosity Burgers experiment solves

ut\+u​ux=ν​ux​x,x∈\[0,1\),t∈\[0,1\],u\_\{t\}\+uu\_\{x\}=\\nu u\_\{xx\},\\qquad x\\in\[0,1\),\\quad t\\in\[0,1\],with

The data are generated on a periodic one\-dimensional grid\. The spatial discretization uses wrap\-around neighbors, corresponding to the periodic identification

u​\(0,t\)=u​\(1,t\),u\(0,t\)=u\(1,t\),with periodic fluxes across the endpoints\.

The input function is the initial condition, sampled as a smooth periodic random Fourier field,

u​\(x,0\)=∑k=16\(ak​sin⁡\(2​π​k​x\)\+bk​cos⁡\(2​π​k​x\)\)\.u\(x,0\)=\\sum\_\{k=1\}^\{6\}\\left\(a\_\{k\}\\sin\(2\\pi kx\)\+b\_\{k\}\\cos\(2\\pi kx\)\\right\)\.The coefficients are sampled with decaying magnitude proportional to1/k1/k, and each initial field is normalized so that its maximum absolute amplitude is bounded\. The reference solutions are generated by a periodic finite\-volume solver using a Rusanov numerical flux for the nonlinear convection term and a periodic centered discretization for the diffusion term, with time step5×10−45\\times 10^\{\-4\}\. This setting is a fixed\-parameter operator\-learning problem:

u​\(x,0\)↦u​\(x,t;ν=0\.01\)\.u\(x,0\)\\mapsto u\(x,t;\\nu=0\.01\)\.
Results\.See Table[10](https://arxiv.org/html/2607.28762#A3.T10)and Figure[21](https://arxiv.org/html/2607.28762#A3.F21)\.

![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/burgers1d_deeponet_best_trainloss_test0.png)\(a\)DeepONet Burgers 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/op_barchart/burgers1d_nfm_deeponet_best_trainloss_test0.png)\(b\)FM\-Operator Burgers 1D

Figure 21:Comparison between DeepONet and FM\-Operator on the fixed\-viscosity periodic Burgers equation\.
##### Parametric Burgers Equation

The parametric Burgers experiment extends the fixed\-viscosity setting by varying both the viscosity and the boundary values:

ut\+u​ux=ν​ux​x,x∈\[−1,1\],t∈\[0,1\],u\_\{t\}\+uu\_\{x\}=\\nu u\_\{xx\},\\qquad x\\in\[\-1,1\],\\quad t\\in\[0,1\],with

ν∼𝒰​\(0\.005,0\.02\),bL,bR∼𝒰​\(−0\.3,0\.3\)\.\\nu\\sim\\mathcal\{U\}\(0\.005,0\.02\),\\qquad b\_\{L\},b\_\{R\}\\sim\\mathcal\{U\}\(\-0\.3,0\.3\)\.The boundary conditions are

u​\(−1,t\)=bL,u​\(1,t\)=bR\.u\(\-1,t\)=b\_\{L\},\\qquad u\(1,t\)=b\_\{R\}\.The initial condition is constructed as a boundary lift plus a zero\-boundary perturbation:

u​\(x,0\)=\(1−z\)​bL\+z​bR\+∑k=16ak​sin⁡\(k​π​z\),z=x\+12\.u\(x,0\)=\(1\-z\)b\_\{L\}\+zb\_\{R\}\+\\sum\_\{k=1\}^\{6\}a\_\{k\}\\sin\(k\\pi z\),\\qquad z=\\frac\{x\+1\}\{2\}\.The perturbation coefficients are sampled with1/k1/kdecay and scaled so that the perturbation amplitude is at most0\.80\.8\. The solution data are generated with the same finite\-difference Burgers solver, while enforcingbLb\_\{L\}andbRb\_\{R\}at the two boundaries during time integration\.

The branch input contains6464initial\-condition sensors together with the normalized viscosity and the two boundary values:

sparam​\(v\)=\[u0​\(x1\),…,u0​\(x64\),ν^,bL,bR\]∈ℝ67\.s\_\{\\mathrm\{param\}\}\(v\)=\\left\[u\_\{0\}\(x\_\{1\}\),\\ldots,u\_\{0\}\(x\_\{64\}\),\\widehat\{\\nu\},b\_\{L\},b\_\{R\}\\right\]\\in\\mathbb\{R\}^\{67\}\.Thus, the model learns the parametric solution operator

𝒢:\(u0,ν,bL,bR\)↦u,𝒢​\(u0,ν,bL,bR\)​\(x,t\)=u​\(x,t\),\\mathcal\{G\}:\(u\_\{0\},\\nu,b\_\{L\},b\_\{R\}\)\\mapsto u,\\qquad\\mathcal\{G\}\(u\_\{0\},\\nu,b\_\{L\},b\_\{R\}\)\(x,t\)=u\(x,t\),
Results\.See Table[11](https://arxiv.org/html/2607.28762#A3.T11)and Figure[22](https://arxiv.org/html/2607.28762#A3.F22)\.

Table 11:Best\-train\-loss test MSE and global relativeL2L\_\{2\}for the Parametric Burgers equation\.ModelParamsBest stepTest MSEGlobal relativeL2L\_\{2\}DeepONet108,161272171\.003×10−31\.003\\times 10^\{\-3\}7\.939×10−27\.939\\times 10^\{\-2\}FM\-DeepONet \(XDeepFM\)453,785279535\.063×10−45\.063\\times 10^\{\-4\}5\.641×10−25\.641\\times 10^\{\-2\}FM\-DeepONet \(NFM\)207,745206051\.503×10−41\.503\\times 10^\{\-4\}3\.073×10−23\.073\\times 10^\{\-2\}FM\-Operator \(XDeepFM\)340,249259965\.639×10−45\.639\\times 10^\{\-4\}5\.953×10−25\.953\\times 10^\{\-2\}FM\-Operator \(NFM\)42,826287541\.718×10−41\.718\\times 10^\{\-4\}3\.286×10−23\.286\\times 10^\{\-2\}

![Refer to caption](https://arxiv.org/html/2607.28762v1/burgers_plots/p_BU/burgers1d_parametric_bcnu_deeponet_best_trainloss_test0.png)\(a\)DeepONet Burgers 1D
![Refer to caption](https://arxiv.org/html/2607.28762v1/burgers_plots/p_BU/burgers1d_parametric_bcnu_nfm_deeponet_best_trainloss_test0.png)\(b\)FM\-Operator Burgers 1D

Figure 22:Comparison between DeepONet and FM\-Operator on Parametric Burgers 1D\.

#### Observations: Conservation Laws and Shocks

Table[10](https://arxiv.org/html/2607.28762#A3.T10)and Figure[17](https://arxiv.org/html/2607.28762#A3.F17)summarize the results\. Two patterns stand out\.

FM\-DeepONet \(NFM\) leads across nonlinear conservation laws\.On cubic conservation, Buckley–Leverett, and periodic Burgers, FM\-DeepONet \(NFM\) achieves the lowest error, reducing the globalL2L\_\{2\}by 39\.06%, 29\.50%, and 59\.68% relative to DeepONet\. The branch–trunk structure augmented with NFM interactions consistently outperforms both the standard DeepONet and the XDeepFM variants\.

FM\-Operator \(NFM\) is competitive with far fewer parameters\.Across the same benchmarks, FM\-Operator \(NFM\) reduces theL2L\_\{2\}error by 33\.03%, 20\.82%, and 34\.20% while using only 42,439–42,826 parameters—roughly 40% of DeepONet’s count \(Table[9](https://arxiv.org/html/2607.28762#A3.T9)\)\. On the parametric Burgers equation, where viscosity and boundary values are jointly varied, FM\-Operator \(NFM\) achieves anL2L\_\{2\}error of3\.286×10−23\.286\\times 10^\{\-2\}with 42,826 parameters, compared to7\.939×10−27\.939\\times 10^\{\-2\}for DeepONet at 108,161 parameters: a 58\.61% error reduction at 39\.60% of the cost \(Table[11](https://arxiv.org/html/2607.28762#A3.T11)\)\.

XDeepFM underperforms despite higher capacity\.FM\-Operator \(XDeepFM\) uses 339,097 parameters but consistently yields the highest errors across conservation\-law benchmarks, including3\.074×10−13\.074\\times 10^\{\-1\}on cubic conservation \(vs\. DeepONet’s2\.319×10−12\.319\\times 10^\{\-1\}\)\. Deeper explicit interaction structures do not automatically translate to better accuracy for these problems\.

Discontinuous transport\.The square\-pulse advection isolates pure linear transport of a discontinuity\. FM\-Operator \(NFM\) achieves an 11\.94%L2L\_\{2\}reduction over DeepONet \(2\.117×10−12\.117\\times 10^\{\-1\}vs\.2\.404×10−12\.404\\times 10^\{\-1\}\)\. Although the Taylor interpretation of Section 4 does not formally apply to discontinuous solutions, the interaction between sensor values, query coordinates, and statistical features appears to help the model track the moving fronts\.

When interactions do not help\.On the Kuramoto–Sivashinsky equation, DeepONet outperforms all FM variants \(globalL2L\_\{2\}of1\.583×10−21\.583\\times 10^\{\-2\}vs\.2\.065×10−22\.065\\times 10^\{\-2\}for the best FM model\)\. The periodic, quasi\-linear dynamics appear not to benefit from explicit cross\-field interactions\. On the LWR traffic flow equation, all models perform comparably, with FM\-DeepONet \(XDeepFM\) marginally ahead\.

#### Observations: Smooth PDEs

TableLABEL:tab:all\_fmo\_split\_nonsplit\_deeponet\-appreports five\-seed statistics across smooth PDE benchmarks spanning 1D to 10D\. Several trends emerge\.

Dimension\-dependent advantage\.FM\-Operator \(NFM\-Split\) outperforms DeepONet on 4D and 5D convection–diffusion, linear advection, and reaction–diffusion, with relativeL2L\_\{2\}reductions of 30–50%\. In 10D, the edge narrows: on Poisson 10D, FM\-Operator \(NFM\-Non\-Split\) and DeepONet are comparable \(9\.489×10−19\.489\\times 10^\{\-1\}vs\.9\.824×10−19\.824\\times 10^\{\-1\}\)\. On Heat 10D and Wave 10D, DeepONet retains a clear advantage\.

Helmholtz and low\-dimensional smooth PDEs favor DeepONet\.On the Helmholtz equation \(1D–5D\) and on 1D heat, Poisson, wave, and linear advection, DeepONet consistently achieves the lowest error\. These problems lack substantial cross\-coordinate dependence, so the interaction mechanism provides no structural benefit, and the added parameters may introduce optimization overhead\.

Split vs\. non\-split input\.FM\-Operator \(NFM\-Split\) and FM\-Operator \(NFM\-Non\-Split\) perform comparably across most benchmarks, with neither variant showing a systematic advantage\. This suggests that the interaction mechanism, rather than the specific field\-partitioning strategy, is the primary contributor to the observed performance differences\.

#### Summary

The results support a conditional view of the interaction mechanism\. When the target mapping involves substantial cross\-field dependence—as in multidimensional smooth PDEs and nonlinear conservation laws where sensor values and query coordinates jointly determine the solution—the NFM bi\-interaction layer provides an effective inductive bias, often at a fraction of the parameter count of a standard DeepONet\. When such dependence is weak \(low\-dimensional smooth problems, Helmholtz, Kuramoto–Sivashinsky\), the mechanism offers no consistent advantage\. The XDeepFM comparison further indicates that deeper interaction hierarchies do not automatically improve upon the simpler NFM\-style pairwise interactions for the benchmarks studied here\.

#### Five\-Seed Statistics

Table 12:Five\-seed global relativeL2L\_\{2\}errors \(mean±\\pmsample standard deviation\)\. The best mean in each row is shown in blue\.EquationFM\-Operator\(NFM\-Split\)FM\-Operator\(NFM\-Non\-Split\)DeepONetNonlinear, shock, and BurgersBuckley–Leverett 1D\(7\.531×10−2\)±\(2\.981×10−3\)\\color\[rgb\]\{0,0,1\}\{\(7\.531\\times 10^\{\-2\}\)\\pm\(2\.981\\times 10^\{\-3\}\)\}\(7\.622×10−2\)±\(5\.044×10−3\)\(7\.622\\times 10^\{\-2\}\)\\pm\(5\.044\\times 10^\{\-3\}\)\(9\.104×10−2\)±\(2\.861×10−3\)\(9\.104\\times 10^\{\-2\}\)\\pm\(2\.861\\times 10^\{\-3\}\)Buckley–Leverett 2D diagonal lift\(7\.310×10−2\)±\(2\.688×10−3\)\\color\[rgb\]\{0,0,1\}\{\(7\.310\\times 10^\{\-2\}\)\\pm\(2\.688\\times 10^\{\-3\}\)\}\(7\.356×10−2\)±\(5\.031×10−3\)\(7\.356\\times 10^\{\-2\}\)\\pm\(5\.031\\times 10^\{\-3\}\)\(8\.920×10−2\)±\(7\.877×10−4\)\(8\.920\\times 10^\{\-2\}\)\\pm\(7\.877\\times 10^\{\-4\}\)Parametric Burgers 1D\(3\.896×10−2\)±\(5\.403×10−3\)\(3\.896\\times 10^\{\-2\}\)\\pm\(5\.403\\times 10^\{\-3\}\)\(3\.894×10−2\)±\(6\.682×10−3\)\\color\[rgb\]\{0,0,1\}\{\(3\.894\\times 10^\{\-2\}\)\\pm\(6\.682\\times 10^\{\-3\}\)\}\(8\.534×10−2\)±\(8\.611×10−3\)\(8\.534\\times 10^\{\-2\}\)\\pm\(8\.611\\times 10^\{\-3\}\)Cubic conservation 1D\(1\.450×10−1\)±\(1\.944×10−2\)\\color\[rgb\]\{0,0,1\}\{\(1\.450\\times 10^\{\-1\}\)\\pm\(1\.944\\times 10^\{\-2\}\)\}\(1\.559×10−1\)±\(1\.302×10−2\)\(1\.559\\times 10^\{\-1\}\)\\pm\(1\.302\\times 10^\{\-2\}\)\(2\.341×10−1\)±\(1\.616×10−2\)\(2\.341\\times 10^\{\-1\}\)\\pm\(1\.616\\times 10^\{\-2\}\)Cubic conservation 2D diagonal lift\(1\.244×10−1\)±\(1\.260×10−2\)\\color\[rgb\]\{0,0,1\}\{\(1\.244\\times 10^\{\-1\}\)\\pm\(1\.260\\times 10^\{\-2\}\)\}\(1\.462×10−1\)±\(1\.197×10−2\)\(1\.462\\times 10^\{\-1\}\)\\pm\(1\.197\\times 10^\{\-2\}\)\(2\.174×10−1\)±\(1\.999×10−2\)\(2\.174\\times 10^\{\-1\}\)\\pm\(1\.999\\times 10^\{\-2\}\)Kuramoto–Sivashinsky 1D\(2\.834×10−2\)±\(1\.816×10−3\)\(2\.834\\times 10^\{\-2\}\)\\pm\(1\.816\\times 10^\{\-3\}\)\(2\.762×10−2\)±\(3\.950×10−3\)\(2\.762\\times 10^\{\-2\}\)\\pm\(3\.950\\times 10^\{\-3\}\)\(1\.830×10−2\)±\(1\.627×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.830\\times 10^\{\-2\}\)\\pm\(1\.627\\times 10^\{\-3\}\)\}Kuramoto–Sivashinsky 2D diagonal lift\(2\.570×10−2\)±\(1\.932×10−3\)\(2\.570\\times 10^\{\-2\}\)\\pm\(1\.932\\times 10^\{\-3\}\)\(2\.501×10−2\)±\(1\.727×10−3\)\(2\.501\\times 10^\{\-2\}\)\\pm\(1\.727\\times 10^\{\-3\}\)\(1\.630×10−2\)±\(1\.536×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.630\\times 10^\{\-2\}\)\\pm\(1\.536\\times 10^\{\-3\}\)\}LWR 1D\(5\.793×10−2\)±\(5\.200×10−3\)\\color\[rgb\]\{0,0,1\}\{\(5\.793\\times 10^\{\-2\}\)\\pm\(5\.200\\times 10^\{\-3\}\)\}\(5\.881×10−2\)±\(4\.556×10−3\)\(5\.881\\times 10^\{\-2\}\)\\pm\(4\.556\\times 10^\{\-3\}\)\(8\.126×10−2\)±\(3\.710×10−3\)\(8\.126\\times 10^\{\-2\}\)\\pm\(3\.710\\times 10^\{\-3\}\)LWR 2D diagonal lift\(6\.417×10−2\)±\(7\.229×10−3\)\(6\.417\\times 10^\{\-2\}\)\\pm\(7\.229\\times 10^\{\-3\}\)\(6\.293×10−2\)±\(6\.013×10−3\)\\color\[rgb\]\{0,0,1\}\{\(6\.293\\times 10^\{\-2\}\)\\pm\(6\.013\\times 10^\{\-3\}\)\}\(8\.798×10−2\)±\(5\.714×10−3\)\(8\.798\\times 10^\{\-2\}\)\\pm\(5\.714\\times 10^\{\-3\}\)Periodic Burgers 1D \(fixed viscosity\)\(8\.548×10−2\)±\(5\.901×10−3\)\\color\[rgb\]\{0,0,1\}\{\(8\.548\\times 10^\{\-2\}\)\\pm\(5\.901\\times 10^\{\-3\}\)\}\(9\.464×10−2\)±\(7\.872×10−3\)\(9\.464\\times 10^\{\-2\}\)\\pm\(7\.872\\times 10^\{\-3\}\)\(1\.356×10−1\)±\(8\.765×10−3\)\(1\.356\\times 10^\{\-1\}\)\\pm\(8\.765\\times 10^\{\-3\}\)Square\-pulse advection 1D\(2\.081×10−1\)±\(5\.286×10−3\)\\color\[rgb\]\{0,0,1\}\{\(2\.081\\times 10^\{\-1\}\)\\pm\(5\.286\\times 10^\{\-3\}\)\}\(2\.135×10−1\)±\(1\.751×10−3\)\(2\.135\\times 10^\{\-1\}\)\\pm\(1\.751\\times 10^\{\-3\}\)\(2\.406×10−1\)±\(4\.072×10−3\)\(2\.406\\times 10^\{\-1\}\)\\pm\(4\.072\\times 10^\{\-3\}\)Square\-pulse advection 2D diagonal lift\(1\.702×10−1\)±\(6\.692×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.702\\times 10^\{\-1\}\)\\pm\(6\.692\\times 10^\{\-3\}\)\}\(1\.708×10−1\)±\(9\.162×10−3\)\(1\.708\\times 10^\{\-1\}\)\\pm\(9\.162\\times 10^\{\-3\}\)\(1\.989×10−1\)±\(6\.298×10−3\)\(1\.989\\times 10^\{\-1\}\)\\pm\(6\.298\\times 10^\{\-3\}\)SmoothConvection–diffusion 1D\(4\.578×10−3\)±\(3\.185×10−3\)\(4\.578\\times 10^\{\-3\}\)\\pm\(3\.185\\times 10^\{\-3\}\)\(4\.018×10−3\)±\(2\.299×10−3\)\(4\.018\\times 10^\{\-3\}\)\\pm\(2\.299\\times 10^\{\-3\}\)\(2\.625×10−3\)±\(1\.373×10−3\)\\color\[rgb\]\{0,0,1\}\{\(2\.625\\times 10^\{\-3\}\)\\pm\(1\.373\\times 10^\{\-3\}\)\}Convection–diffusion 2D\(7\.265×10−3\)±\(4\.449×10−3\)\(7\.265\\times 10^\{\-3\}\)\\pm\(4\.449\\times 10^\{\-3\}\)\(6\.496×10−3\)±\(2\.763×10−3\)\\color\[rgb\]\{0,0,1\}\{\(6\.496\\times 10^\{\-3\}\)\\pm\(2\.763\\times 10^\{\-3\}\)\}\(9\.717×10−3\)±\(7\.468×10−3\)\(9\.717\\times 10^\{\-3\}\)\\pm\(7\.468\\times 10^\{\-3\}\)Convection–diffusion 4D\(9\.925×10−3\)±\(2\.108×10−3\)\\color\[rgb\]\{0,0,1\}\{\(9\.925\\times 10^\{\-3\}\)\\pm\(2\.108\\times 10^\{\-3\}\)\}\(1\.163×10−2\)±\(3\.150×10−3\)\(1\.163\\times 10^\{\-2\}\)\\pm\(3\.150\\times 10^\{\-3\}\)\(1\.874×10−2\)±\(1\.728×10−3\)\(1\.874\\times 10^\{\-2\}\)\\pm\(1\.728\\times 10^\{\-3\}\)Convection–diffusion 5D\(1\.444×10−2\)±\(2\.958×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.444\\times 10^\{\-2\}\)\\pm\(2\.958\\times 10^\{\-3\}\)\}\(1\.462×10−2\)±\(2\.629×10−3\)\(1\.462\\times 10^\{\-2\}\)\\pm\(2\.629\\times 10^\{\-3\}\)\(3\.095×10−2\)±\(3\.196×10−3\)\(3\.095\\times 10^\{\-2\}\)\\pm\(3\.196\\times 10^\{\-3\}\)Heat 1D\(8\.130×10−3\)±\(8\.362×10−4\)\(8\.130\\times 10^\{\-3\}\)\\pm\(8\.362\\times 10^\{\-4\}\)\(8\.670×10−3\)±\(2\.327×10−3\)\(8\.670\\times 10^\{\-3\}\)\\pm\(2\.327\\times 10^\{\-3\}\)\(4\.430×10−3\)±\(4\.244×10−4\)\\color\[rgb\]\{0,0,1\}\{\(4\.430\\times 10^\{\-3\}\)\\pm\(4\.244\\times 10^\{\-4\}\)\}Heat 4D\(6\.594×10−2\)±\(3\.224×10−3\)\(6\.594\\times 10^\{\-2\}\)\\pm\(3\.224\\times 10^\{\-3\}\)\(5\.886×10−2\)±\(2\.610×10−3\)\\color\[rgb\]\{0,0,1\}\{\(5\.886\\times 10^\{\-2\}\)\\pm\(2\.610\\times 10^\{\-3\}\)\}\(6\.261×10−2\)±\(3\.406×10−3\)\(6\.261\\times 10^\{\-2\}\)\\pm\(3\.406\\times 10^\{\-3\}\)Heat 10D\(2\.380×10−1\)±\(9\.563×10−3\)\(2\.380\\times 10^\{\-1\}\)\\pm\(9\.563\\times 10^\{\-3\}\)\(2\.330×10−1\)±\(1\.970×10−2\)\(2\.330\\times 10^\{\-1\}\)\\pm\(1\.970\\times 10^\{\-2\}\)\(1\.877×10−1\)±\(2\.053×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.877\\times 10^\{\-1\}\)\\pm\(2\.053\\times 10^\{\-3\}\)\}Helmholtz 1D\(5\.520×10−2\)±\(5\.062×10−2\)\(5\.520\\times 10^\{\-2\}\)\\pm\(5\.062\\times 10^\{\-2\}\)\(6\.186×10−2\)±\(2\.686×10−2\)\(6\.186\\times 10^\{\-2\}\)\\pm\(2\.686\\times 10^\{\-2\}\)\(2\.285×10−2\)±\(1\.686×10−2\)\\color\[rgb\]\{0,0,1\}\{\(2\.285\\times 10^\{\-2\}\)\\pm\(1\.686\\times 10^\{\-2\}\)\}Helmholtz 2D\(5\.246×10−2\)±\(1\.894×10−2\)\(5\.246\\times 10^\{\-2\}\)\\pm\(1\.894\\times 10^\{\-2\}\)\(7\.258×10−2\)±\(3\.372×10−2\)\(7\.258\\times 10^\{\-2\}\)\\pm\(3\.372\\times 10^\{\-2\}\)\(7\.891×10−3\)±\(3\.863×10−3\)\\color\[rgb\]\{0,0,1\}\{\(7\.891\\times 10^\{\-3\}\)\\pm\(3\.863\\times 10^\{\-3\}\)\}Helmholtz 4D\(1\.492×10−1\)±\(7\.721×10−2\)\(1\.492\\times 10^\{\-1\}\)\\pm\(7\.721\\times 10^\{\-2\}\)\(2\.057×10−1\)±\(1\.044×10−1\)\(2\.057\\times 10^\{\-1\}\)\\pm\(1\.044\\times 10^\{\-1\}\)\(2\.746×10−2\)±\(1\.229×10−2\)\\color\[rgb\]\{0,0,1\}\{\(2\.746\\times 10^\{\-2\}\)\\pm\(1\.229\\times 10^\{\-2\}\)\}Helmholtz 5D\(3\.043×10−1\)±\(1\.619×10−1\)\(3\.043\\times 10^\{\-1\}\)\\pm\(1\.619\\times 10^\{\-1\}\)\(6\.190×10−1\)±\(4\.890×10−1\)\(6\.190\\times 10^\{\-1\}\)\\pm\(4\.890\\times 10^\{\-1\}\)\(8\.894×10−2\)±\(7\.995×10−2\)\\color\[rgb\]\{0,0,1\}\{\(8\.894\\times 10^\{\-2\}\)\\pm\(7\.995\\times 10^\{\-2\}\)\}Linear advection 1D\(5\.520×10−3\)±\(1\.675×10−3\)\(5\.520\\times 10^\{\-3\}\)\\pm\(1\.675\\times 10^\{\-3\}\)\(5\.210×10−3\)±\(2\.327×10−3\)\(5\.210\\times 10^\{\-3\}\)\\pm\(2\.327\\times 10^\{\-3\}\)\(3\.370×10−3\)±\(1\.439×10−3\)\\color\[rgb\]\{0,0,1\}\{\(3\.370\\times 10^\{\-3\}\)\\pm\(1\.439\\times 10^\{\-3\}\)\}Linear advection 2D\(4\.548×10−3\)±\(5\.949×10−4\)\\color\[rgb\]\{0,0,1\}\{\(4\.548\\times 10^\{\-3\}\)\\pm\(5\.949\\times 10^\{\-4\}\)\}\(4\.813×10−3\)±\(9\.255×10−4\)\(4\.813\\times 10^\{\-3\}\)\\pm\(9\.255\\times 10^\{\-4\}\)\(4\.585×10−3\)±\(5\.461×10−4\)\(4\.585\\times 10^\{\-3\}\)\\pm\(5\.461\\times 10^\{\-4\}\)Linear advection 4D\(1\.061×10−2\)±\(2\.714×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.061\\times 10^\{\-2\}\)\\pm\(2\.714\\times 10^\{\-3\}\)\}\(1\.123×10−2\)±\(2\.339×10−3\)\(1\.123\\times 10^\{\-2\}\)\\pm\(2\.339\\times 10^\{\-3\}\)\(1\.965×10−2\)±\(1\.716×10−3\)\(1\.965\\times 10^\{\-2\}\)\\pm\(1\.716\\times 10^\{\-3\}\)Linear advection 5D\(1\.251×10−2\)±\(1\.023×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.251\\times 10^\{\-2\}\)\\pm\(1\.023\\times 10^\{\-3\}\)\}\(1\.268×10−2\)±\(7\.592×10−4\)\(1\.268\\times 10^\{\-2\}\)\\pm\(7\.592\\times 10^\{\-4\}\)\(3\.043×10−2\)±\(2\.951×10−3\)\(3\.043\\times 10^\{\-2\}\)\\pm\(2\.951\\times 10^\{\-3\}\)Poisson 1D\(9\.839×10−3\)±\(2\.546×10−3\)\(9\.839\\times 10^\{\-3\}\)\\pm\(2\.546\\times 10^\{\-3\}\)\(9\.855×10−3\)±\(1\.354×10−3\)\(9\.855\\times 10^\{\-3\}\)\\pm\(1\.354\\times 10^\{\-3\}\)\(2\.999×10−3\)±\(1\.338×10−3\)\\color\[rgb\]\{0,0,1\}\{\(2\.999\\times 10^\{\-3\}\)\\pm\(1\.338\\times 10^\{\-3\}\)\}Poisson 4D\(2\.544×10−1\)±\(2\.339×10−2\)\(2\.544\\times 10^\{\-1\}\)\\pm\(2\.339\\times 10^\{\-2\}\)\(3\.762×10−1\)±\(4\.894×10−2\)\(3\.762\\times 10^\{\-1\}\)\\pm\(4\.894\\times 10^\{\-2\}\)\(6\.796×10−2\)±\(2\.381×10−3\)\\color\[rgb\]\{0,0,1\}\{\(6\.796\\times 10^\{\-2\}\)\\pm\(2\.381\\times 10^\{\-3\}\)\}Poisson 10D\(1\.051×100\)±\(2\.703×10−2\)\(1\.051\\times 10^\{0\}\)\\pm\(2\.703\\times 10^\{\-2\}\)\(9\.489×10−1\)±\(9\.658×10−3\)\\color\[rgb\]\{0,0,1\}\{\(9\.489\\times 10^\{\-1\}\)\\pm\(9\.658\\times 10^\{\-3\}\)\}\(9\.824×10−1\)±\(6\.985×10−3\)\(9\.824\\times 10^\{\-1\}\)\\pm\(6\.985\\times 10^\{\-3\}\)Reaction–diffusion 1D\(3\.391×10−2\)±\(4\.208×10−2\)\\color\[rgb\]\{0,0,1\}\{\(3\.391\\times 10^\{\-2\}\)\\pm\(4\.208\\times 10^\{\-2\}\)\}\(3\.663×10−2\)±\(5\.754×10−2\)\(3\.663\\times 10^\{\-2\}\)\\pm\(5\.754\\times 10^\{\-2\}\)\(3\.636×10−2\)±\(6\.502×10−2\)\(3\.636\\times 10^\{\-2\}\)\\pm\(6\.502\\times 10^\{\-2\}\)Reaction–diffusion 2D\(8\.916×10−3\)±\(7\.069×10−3\)\(8\.916\\times 10^\{\-3\}\)\\pm\(7\.069\\times 10^\{\-3\}\)\(8\.849×10−3\)±\(6\.364×10−3\)\(8\.849\\times 10^\{\-3\}\)\\pm\(6\.364\\times 10^\{\-3\}\)\(8\.024×10−3\)±\(8\.438×10−3\)\\color\[rgb\]\{0,0,1\}\{\(8\.024\\times 10^\{\-3\}\)\\pm\(8\.438\\times 10^\{\-3\}\)\}Reaction–diffusion 4D\(8\.601×10−3\)±\(1\.057×10−3\)\\color\[rgb\]\{0,0,1\}\{\(8\.601\\times 10^\{\-3\}\)\\pm\(1\.057\\times 10^\{\-3\}\)\}\(9\.855×10−3\)±\(2\.126×10−3\)\(9\.855\\times 10^\{\-3\}\)\\pm\(2\.126\\times 10^\{\-3\}\)\(1\.595×10−2\)±\(1\.743×10−3\)\(1\.595\\times 10^\{\-2\}\)\\pm\(1\.743\\times 10^\{\-3\}\)Reaction–diffusion 5D\(1\.094×10−2\)±\(1\.095×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.094\\times 10^\{\-2\}\)\\pm\(1\.095\\times 10^\{\-3\}\)\}\(1\.235×10−2\)±\(2\.216×10−3\)\(1\.235\\times 10^\{\-2\}\)\\pm\(2\.216\\times 10^\{\-3\}\)\(2\.917×10−2\)±\(4\.017×10−3\)\(2\.917\\times 10^\{\-2\}\)\\pm\(4\.017\\times 10^\{\-3\}\)Wave 1D\(2\.564×10−2\)±\(3\.324×10−3\)\(2\.564\\times 10^\{\-2\}\)\\pm\(3\.324\\times 10^\{\-3\}\)\(2\.671×10−2\)±\(7\.214×10−3\)\(2\.671\\times 10^\{\-2\}\)\\pm\(7\.214\\times 10^\{\-3\}\)\(1\.154×10−2\)±\(1\.510×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.154\\times 10^\{\-2\}\)\\pm\(1\.510\\times 10^\{\-3\}\)\}Wave 4D\(1\.848×10−1\)±\(1\.009×10−2\)\(1\.848\\times 10^\{\-1\}\)\\pm\(1\.009\\times 10^\{\-2\}\)\(1\.831×10−1\)±\(8\.104×10−3\)\(1\.831\\times 10^\{\-1\}\)\\pm\(8\.104\\times 10^\{\-3\}\)\(1\.470×10−1\)±\(7\.517×10−3\)\\color\[rgb\]\{0,0,1\}\{\(1\.470\\times 10^\{\-1\}\)\\pm\(7\.517\\times 10^\{\-3\}\)\}Wave 10D\(3\.803×10−1\)±\(4\.611×10−2\)\\color\[rgb\]\{0,0,1\}\{\(3\.803\\times 10^\{\-1\}\)\\pm\(4\.611\\times 10^\{\-2\}\)\}\(6\.453×10−1\)±\(1\.256×10−1\)\(6\.453\\times 10^\{\-1\}\)\\pm\(1\.256\\times 10^\{\-1\}\)\(4\.930×10−1\)±\(4\.504×10−2\)\(4\.930\\times 10^\{\-1\}\)\\pm\(4\.504\\times 10^\{\-2\}\)Table 13:Number of trainable network parameters for NFM\-Split, NFM\-Non\-Split, and DeepONet on each benchmark\. Parameter counts are identical across the five random seeds\.EquationFM\-Operator\(NFM\-Split\)FM\-Operator\(NFM\-Non\-Split\)DeepONetNonlinear, shock, and BurgersBuckley–Leverett 1D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Buckley–Leverett 2D diagonal lift42,82442\{,\}82442,56842\{,\}568107,905107\{,\}905Parametric Burgers 1D42,95442\{,\}95442,82642\{,\}826108,161108\{,\}161Cubic conservation 1D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Cubic conservation 2D diagonal lift42,82442\{,\}82442,56842\{,\}568107,905107\{,\}905Kuramoto–Sivashinsky 1D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Kuramoto–Sivashinsky 2D diagonal lift42,82442\{,\}82442,56842\{,\}568107,905107\{,\}905LWR 1D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777LWR 2D diagonal lift42,82442\{,\}82442,56842\{,\}568107,905107\{,\}905Periodic Burgers 1D \(fixed viscosity\)42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Square\-pulse advection 1D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Square\-pulse advection 2D diagonal lift42,82442\{,\}82442,56842\{,\}568107,905107\{,\}905SmoothConvection–diffusion 1D42,31042\{,\}31042,31042\{,\}310107,649107\{,\}649Convection–diffusion 2D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Convection–diffusion 4D43,08143\{,\}08142,69742\{,\}697108,033108\{,\}033Convection–diffusion 5D43,33843\{,\}33842,82642\{,\}826108,161108\{,\}161Heat 1D42,56742\{,\}56742,69642\{,\}696107,777107\{,\}777Heat 4D43,08143\{,\}08142,95442\{,\}954108,033108\{,\}033Heat 10D44,62344\{,\}62343,72843\{,\}728108,801108\{,\}801Helmholtz 1D42,31042\{,\}31042,31042\{,\}310107,649107\{,\}649Helmholtz 2D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Helmholtz 4D43,08143\{,\}08142,69742\{,\}697108,033108\{,\}033Helmholtz 5D43,33843\{,\}33842,82642\{,\}826108,161108\{,\}161Linear advection 1D42,31042\{,\}31042,31042\{,\}310107,649107\{,\}649Linear advection 2D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Linear advection 4D43,08143\{,\}08142,69742\{,\}697108,033108\{,\}033Linear advection 5D43,33843\{,\}33842,82642\{,\}826108,161108\{,\}161Poisson 1D42,31042\{,\}31042,31042\{,\}310107,649107\{,\}649Poisson 4D43,08143\{,\}08142,69742\{,\}697108,033108\{,\}033Poisson 10D44,62344\{,\}62343,47143\{,\}471108,801108\{,\}801Reaction–diffusion 1D42,31042\{,\}31042,31042\{,\}310107,649107\{,\}649Reaction–diffusion 2D42,56742\{,\}56742,43942\{,\}439107,777107\{,\}777Reaction–diffusion 4D43,08143\{,\}08142,69742\{,\}697108,033108\{,\}033Reaction–diffusion 5D43,33843\{,\}33842,82642\{,\}826108,161108\{,\}161Wave 1D42,56742\{,\}56742,69642\{,\}696107,777107\{,\}777Wave 4D43,08143\{,\}08142,95442\{,\}954108,033108\{,\}033Wave 10D44,62344\{,\}62343,72843\{,\}728108,801108\{,\}801TableLABEL:tab:all\_fmo\_split\_nonsplit\_deeponet\-appsummarizes the five\-seed comparison among NFM\-Split, NFM\-Non\-Split, and DeepONet over all 37 benchmarks\. Taking the better mean of the two FM\-Operator variants on each benchmark, FM\-Operator achieves a lower mean relativeL2L\_\{2\}error than DeepONet on 22 benchmarks, while DeepONet performs best on the remaining 15\. More specifically, NFM\-Split gives the lowest mean in 17 rows, NFM\-Non\-Split in five rows, and DeepONet in 15 rows\. Thus, the results show that the effectiveness of FM\-Operator is problem\-dependent rather than uniform across all operator\-learning tasks\.

The clearest advantage of FM\-Operator appears in the nonlinear, shock, and Burgers group, where it achieves the lowest mean error on 10 of the 12 benchmarks\. Relative to DeepONet, the best FM\-Operator variant reduces the mean error by17\.3%17\.3\\%and18\.0%18\.0\\%on the 1D and 2D diagonal\-lift Buckley–Leverett problems, respectively; by38\.1%38\.1\\%and42\.8%42\.8\\%on the corresponding cubic conservation problems; by28\.7%28\.7\\%and28\.5%28\.5\\%on LWR; and by13\.5%13\.5\\%and14\.4%14\.4\\%on square\-pulse advection\. It also yields reductions of37\.0%37\.0\\%on periodic Burgers and54\.4%54\.4\\%on parametric Burgers\. The two Kuramoto–Sivashinsky benchmarks are the exceptions: the best FM\-Operator mean is50\.9%50\.9\\%higher than DeepONet in 1D and53\.4%53\.4\\%higher for the 2D diagonal lift\. These results support FM\-Operator as a strong alternative for several nonlinear and discontinuity\-dominated operators, but not as a uniformly superior method\.

For the smooth benchmark families, the comparison is more mixed and often dimension\-dependent\. In convection–diffusion, FM\-Operator is53\.1%53\.1\\%worse than DeepONet in 1D, but lowers the mean error by33\.1%33\.1\\%,47\.0%47\.0\\%, and53\.3%53\.3\\%in 2D, 4D, and 5D, respectively\. Linear advection shows a similar pattern: FM\-Operator is54\.6%54\.6\\%worse in 1D, is only marginally better in 2D \(0\.8%0\.8\\%\), and improves the 4D and 5D results by46\.0%46\.0\\%and58\.9%58\.9\\%\. For reaction–diffusion, FM\-Operator improves the 1D, 4D, and 5D means by6\.7%6\.7\\%,46\.1%46\.1\\%, and62\.5%62\.5\\%, respectively, but is10\.3%10\.3\\%worse in 2D\. These families suggest that FM\-Operator can become more effective as the dimensionality and coordinate interactions increase, although the trend is not universal\.

Indeed, FM\-Operator wins only the 4D heat benchmark, only the 10D Poisson benchmark, and only the 10D wave benchmark within those families, with improvements of6\.0%6\.0\\%,3\.4%3\.4\\%, and22\.9%22\.9\\%, respectively\. DeepONet achieves the lowest mean on all four Helmholtz benchmarks and is substantially more accurate in each case\. Overall, the results indicate that FM\-Operator provides a specialized inductive bias whose benefit depends on the structure of the target operator\. NFM\-Split is the stronger variant more often, but NFM\-Non\-Split gives the best result on five benchmarks, so neither formulation dominates in every setting\.

Taken together, the results provide evidence that FM\-Operator is particularly well suited to certain high\-dimensional operator\-learning problems\. Across several matched equation families, including convection–diffusion, linear advection, reaction–diffusion, Poisson, and wave equations, its relative performance generally improves as the input dimension increases\. This pattern is consistent with the hypothesis that explicitly modeling feature interactions becomes more beneficial when the coordinate structure is higher\-dimensional and more strongly coupled\. However, the non\-monotone heat results and the consistently unfavorable Helmholtz results show that dimensionality alone does not determine performance\. The advantage therefore appears to depend jointly on dimension and operator structure\.

#### Comparison with Shift\-DeepONet

Shift\-DeepONet\(Hadorn[2022](https://arxiv.org/html/2607.28762#bib.bib22)\)extends the standard DeepONet by allowing its trunk basis functions to adapt to the input function\. A conventional DeepONet represents the learned operator as

𝒢^​\(s\)​\(y\)=∑k=1pbk​\(s\)​tk​\(y\),\\widehat\{\\mathcal\{G\}\}\(s\)\(y\)=\\sum\_\{k=1\}^\{p\}b\_\{k\}\(s\)\\,t\_\{k\}\(y\),\(12\)wheressdenotes the sampled input function and the trunk basistk​\(y\)t\_\{k\}\(y\)is shared across all inputs\. Shift\-DeepONet instead evaluates each basis function at an input\-dependent affine transformation of the query coordinate:

𝒢^shift​\(s\)​\(y\)=∑k=1pbk​\(s\)​tk​\(Ak​\(s\)​y\+γk​\(s\)\)\.\\widehat\{\\mathcal\{G\}\}\_\{\\mathrm\{shift\}\}\(s\)\(y\)=\\sum\_\{k=1\}^\{p\}b\_\{k\}\(s\)\\,t\_\{k\}\\\!\\left\(A\_\{k\}\(s\)y\+\\gamma\_\{k\}\(s\)\\right\)\.\(13\)The transformationsAk​\(s\)A\_\{k\}\(s\)and shiftsγk​\(s\)\\gamma\_\{k\}\(s\)allow the learned basis functions to move and deform with the input\. This is particularly useful when sharp fronts or discontinuities occur at different locations, since a fixed trunk basis would otherwise have to represent all possible front positions simultaneously\. Hadorn reports substantial improvements over standard DeepONet on problems with moving discontinuities and sharp transitions\.

We compare FM\-Operator \(NFM\) with this affine\-coordinate formulation of Shift\-DeepONet on the viscous Burgers equation withν=0\.001\\nu=0\.001,x∈\[−1,1\]x\\in\[\-1,1\],t∈\[0,1\]t\\in\[0,1\], and homogeneous Dirichlet boundary conditions\. Both models use the same three random seeds, training and test functions, WENO5 reference solutions, optimization budget, batch construction, and checkpoint\-selection rule\. Table[14](https://arxiv.org/html/2607.28762#A3.T14)reports the mean and sample standard deviation of the global relativeL2L\_\{2\}error\.

Table 14:FM\-Operator \(NFM\) and Shift\-DeepONet on the viscous Burgers equation withν=0\.001\\nu=0\.001\. Results are reported over three seeds as mean±\\pmsample standard deviation\.ModelParametersGlobal relativeL2L\_\{2\}Error reductionShift\-DeepONet∼\\sim141K0\.1132±0\.01570\.1132\\pm 0\.0157—FM\-Operator \(NFM\)∼\\sim43K0\.0815±0\.00240\.0815\\pm 0\.002428\.0%FM\-Operator reduces the global relativeL2L\_\{2\}error by 28\.0% relative to Shift\-DeepONet while using approximately 70% fewer parameters \(∼\\sim43K versus∼\\sim141K\)\. Shift\-DeepONet devotes additional capacity to input\-dependent transformations of its trunk basis\. FM\-Operator instead directly models pairwise interactions among the sampled input function, the spatial and temporal coordinates, and function\-level statistics\. On this low\-viscosity Burgers problem, these cross\-field interactions provide both better accuracy and substantially greater parameter efficiency\.

#### Parameter\-Controlled Ablations

FM\-Operator uses 42,567 parameters in these benchmarks, compared with 107,777 for the standard DeepONet\. This sizeable gap raises a natural question: does FM\-Operator benefit from its interaction structure, or simply from having fewer parameters? We examine this question from both directions\. First, we reduce the size of DeepONet and compare it with FM\-Operator\. We then compare the standard DeepONet with a compact FM\-DeepONet of similar parameter count\. Unless stated otherwise, all reported values are global relativeL2L\_\{2\}errors from the checkpoint with the lowest training loss\.

Short\-DeepONet\.We reduce the branch and trunk widths of DeepONet, yielding a 58,241\-parameter model\. This substantially narrows the capacity gap to FM\-Operator \(NFM\), which has 42,567 parameters\. Table[15](https://arxiv.org/html/2607.28762#A3.T15)reports the mean over three independent seeds on five nonlinear benchmarks\.

Table 15:FM\-Operator \(NFM\) vs\. reduced\-capacity DeepONet\. Results are three\-seed means of the global relativeL2L\_\{2\}error\.Equationshort\-DeepONetL2L\_\{2\}FM\-OperatorL2L\_\{2\}ReductionCubic conservation0\.23270\.142238\.9%Buckley–Leverett0\.11340\.075533\.5%Square\-pulse advection0\.32490\.203837\.3%Periodic Burgers0\.17150\.078754\.1%Kuramoto–Sivashinsky0\.09000\.027769\.2%FM\-Operator achieves a lower error on every benchmark while using 26\.9% fewer parameters than short\-DeepONet\. The reduction ranges from 33\.5% on Buckley–Leverett to 69\.2% on Kuramoto–Sivashinsky\. Thus, narrowing the parameter gap does not close the accuracy gap: in this regime, the explicit cross\-field interactions provide a more effective representation than adding capacity to a conventional branch–trunk model\.

Compact FM\-DeepONet\.We next compare the standard DeepONet \(107,777 parameters\) with a compact FM\-DeepONet \(NFM, 124,929 parameters\), a difference of 15\.9%\. The DeepONet results are the five\-seed means from the main study, while the compact FM\-DeepONet results are three\-seed means from this ablation\.

Table 16:Compact FM\-DeepONet \(NFM\) vs\. standard DeepONet at similar parameter counts\. Entries are mean global relativeL2L\_\{2\}errors\.EquationDeepONetL2L\_\{2\}FM\-DeepONetL2L\_\{2\}ReductionCubic conservation0\.23410\.152135\.0%Buckley–Leverett0\.09100\.068824\.4%Square\-pulse advection0\.24060\.210212\.6%Periodic Burgers0\.13560\.078042\.5%Kuramoto–Sivashinsky0\.01830\.0203−\-10\.7%FM\-DeepONet outperforms standard DeepONet on four of five benchmarks, with the largest reductions on periodic Burgers \(42\.5%\) and cubic conservation \(35\.0%\)\. Kuramoto–Sivashinsky is the only exception, with FM\-DeepONet producing a 10\.7% higher error\. This exception is important: the interaction module is helpful on most of these nonlinear benchmarks, but its benefit is not universal, and the standard branch–trunk factorization remains better suited to the Kuramoto–Sivashinsky dynamics in this setting\.

DeepONet with statistics\.FM\-Operator and FM\-DeepONet receive three additional scalar inputs \(mean, standard deviation, and maximum absolute value of the sensor vector\) that are not supplied to the baseline DeepONet\. We therefore append the same statistics to the DeepONet branch input and leave the rest of the architecture unchanged\. Table[17](https://arxiv.org/html/2607.28762#A3.T17)reports the results on eight benchmarks\.

Table 17:DeepONet with vs\. without sensor statistics\.EquationDeepONetL2L\_\{2\}DeepONet\+\+statsL2L\_\{2\}ChangeBurgers \(Dirichlet\)7\.644×10−27\.644\\times 10^\{\-2\}7\.105×10−27\.105\\times 10^\{\-2\}−\-7\.1%Cubic conservation2\.319×10−12\.319\\times 10^\{\-1\}2\.291×10−12\.291\\times 10^\{\-1\}−\-1\.2%Buckley–Leverett8\.649×10−28\.649\\times 10^\{\-2\}8\.637×10−28\.637\\times 10^\{\-2\}−\-0\.1%Square\-pulse advection2\.404×10−12\.404\\times 10^\{\-1\}2\.470×10−12\.470\\times 10^\{\-1\}\+\+2\.7%Periodic Burgers1\.282×10−11\.282\\times 10^\{\-1\}1\.344×10−11\.344\\times 10^\{\-1\}\+\+4\.8%Parametric Burgers7\.939×10−27\.939\\times 10^\{\-2\}8\.290×10−28\.290\\times 10^\{\-2\}\+\+4\.4%Kuramoto–Sivashinsky1\.583×10−21\.583\\times 10^\{\-2\}1\.807×10−21\.807\\times 10^\{\-2\}\+\+14\.1%LWR traffic flow8\.251×10−28\.251\\times 10^\{\-2\}8\.252×10−28\.252\\times 10^\{\-2\}\+\+0\.0%Adding these statistics does not produce a consistent gain\. The mean of the eight per\-benchmark error ratios is1\.0221\.022\(statistics/baseline\), and five benchmarks become less accurate\. The improvements on the other three benchmarks are all below 8%, far smaller than the reductions observed for the interaction\-based models on the corresponding problems\. The three summary statistics alone therefore cannot explain the performance gains\.

Together, these ablations tell a consistent story\. A smaller DeepONet does not match FM\-Operator, a compact FM\-DeepONet remains stronger than the standard DeepONet on four of the five tests, and simply appending the sensor statistics does not reproduce the improvement\. The gains are therefore best attributed to how the models construct cross\-field interactions, rather than to parameter count or to the three auxiliary statistics alone\.

## Appendix DSensitivity of the SPINN Comparison to Solution Structure

This appendix compares FM\-PINN with Separable PINN \(SPINN\) under two manufactured\-solution families\. The comparison is motivated by the different inductive biases of the two architectures\. SPINN represents the solution as a sum of rank\-one, dimension\-wise separable factors, whereas FM\-PINN explicitly mixes coordinate fields before nonlinear decoding\. Consequently, a benchmark whose exact solution is already a short sum of separable products may favor SPINN, while non\-axis\-aligned coordinate couplings provide a more direct test of the interaction mechanism\.

### Matched Experimental Protocol

We denote the original smooth benchmark by*Setting S*and the additional quadratic\-phase benchmark by*Setting Q*\. The differential operators are identical in both settings\. OnΩ=\[0,1\]d\\Omega=\[0,1\]^\{d\}, we solve

−0\.05​Δ​u\+𝜷⋅∇u=fCD,−0\.05​Δ​u\+\(u3−u\)=fRD,\-0\.05\\Delta u\+\\boldsymbol\{\\beta\}\\cdot\\nabla u=f\_\{\\mathrm\{CD\}\},\\qquad\-0\.05\\Delta u\+\(u^\{3\}\-u\)=f\_\{\\mathrm\{RD\}\},where𝜷=\(1,−0\.7\)\\boldsymbol\{\\beta\}=\(1,\-0\.7\)in two dimensions and𝜷=\(1,−0\.7,0\.5,−0\.3\)\\boldsymbol\{\\beta\}=\(1,\-0\.7,0\.5,\-0\.3\)in four dimensions\. In each case, the forcing term is obtained by substituting the manufactured exact solution into the corresponding operator\.

Setting S uses the smooth solution already defined in Appendix C,

uS∗​\(𝒙\)=∑m=15cm​∏j=1dsin⁡\(km,j​π​xj\)\.u\_\{\\mathrm\{S\}\}^\{\*\}\(\\boldsymbol\{x\}\)=\\sum\_\{m=1\}^\{5\}c\_\{m\}\\prod\_\{j=1\}^\{d\}\\sin\(k\_\{m,j\}\\pi x\_\{j\}\)\.It is a sum of at most five rank\-one separable functions and therefore has CP rank at most five after tensor\-grid discretization\. This structure is closely aligned with the rank\-128 SPINN representation\.

Setting Q retains a zero\-boundary trigonometric envelope but introduces quadratic coordinate phases\. Letzi=xi−1/2z\_\{i\}=x\_\{i\}\-1/2\. In two dimensions,

q1=z1​z2,q2=z12−z222,q\_\{1\}=z\_\{1\}z\_\{2\},\\qquad q\_\{2\}=\\frac\{z\_\{1\}^\{2\}\-z\_\{2\}^\{2\}\}\{2\},and in four dimensions,

q1=z1​z2\+z3​z4,q2=z1​z3\+z2​z4\.q\_\{1\}=z\_\{1\}z\_\{2\}\+z\_\{3\}z\_\{4\},\\qquad q\_\{2\}=z\_\{1\}z\_\{3\}\+z\_\{2\}z\_\{4\}\.The exact solution is

uQ∗​\(𝒙\)=∏i=1dsin⁡\(π​xi\)​\[1\+0\.6​sin⁡\(2​π​q1\)\+0\.4​cos⁡\(2​π​q2\)\]\.u\_\{\\mathrm\{Q\}\}^\{\*\}\(\\boldsymbol\{x\}\)=\\prod\_\{i=1\}^\{d\}\\sin\(\\pi x\_\{i\}\)\\left\[1\+0\.6\\sin\(2\\pi q\_\{1\}\)\+0\.4\\cos\(2\\pi q\_\{2\}\)\\right\]\.\(14\)Unlike Setting S, Eq\. \([14](https://arxiv.org/html/2607.28762#A4.E14)\) contains non\-axis\-aligned coordinate products inside nonlinear phase functions and is not a prescribed finite sum of coordinate\-wise rank\-one terms\.

Table[18](https://arxiv.org/html/2607.28762#A4.T18)summarizes the numerical protocol\. All principal settings are held fixed\. Setting Q additionally aligns the mini\-batch index sequence across models, whereas Setting S uses the same LHS pool with model\-specific deterministic mini\-batch streams\. FM\-PINN uses a 64\-dimensional interaction feature and a width\-128, depth\-4 decoder\. SPINN uses one width\-64, depth\-4 body network per coordinate and CP rank 128\.

Table 18:Protocol for the two SPINN comparison settings\.ConfigurationSetting SSetting QRandom seed9191Training steps30,00030,000LHS interior points10,00010,000Boundary points400400Interior/boundary batch1,024/4001,024/400OptimizerAdamWAdamWLearning rate/weight decay10−3/10−610^\{\-3\}/10^\{\-6\}10−3/10−610^\{\-3\}/10^\{\-6\}Loss weights \(PDE/boundary\)1/11/1Evaluation grid \(2D/4D\)1282/244128^\{2\}/24^\{4\}1282/244128^\{2\}/24^\{4\}SPINN rank/width/depth128/64/4128/64/4Sampling pool across modelsSharedSharedMini\-batch indices across modelsModel\-specificSharedExact\-solution structure≤5\\leq 5separable termsQuadratic phases
### Results

TableLABEL:tab:spinn\_solution\_structure\_resultsreports parameter count, wall\-clock time, and relativeL2L\_\{2\}error\. “Best” denotes the test error of the checkpoint selected by the minimum stochastic training loss; it is not a checkpoint selected on the test grid\. “Final” evaluates the parameters at the fixed 30,000\-step budget\. The Setting\-S SPINN entries are the four records in the original SPINN comparison file; the accompanying FM\-PINN records are the completed runs under the same seed and protocol\. Setting Q is the subsequent 16\-run comparison\. We retain the first\-order FM variant in the table for completeness, but focus the analysis on FM\-PINN without the explicit linear output branch because the first\-order branch gives no consistent gain and can increase runtime\.

Table 19:Accuracy and runtime under separable \(S\) and quadratic\-phase \(Q\) manufactured solutions\. The lowest best\-checkpoint error within each PDE/setting block is shown in bold\.SettingPDEModelParamsBest rel\.L2L\_\{2\}Final rel\.L2L\_\{2\}SCD–2DFM\-PINN58,2411\.031×10−31\.031\\text\{\\times\}\{10\}^\{\-3\}8\.269×10−38\.269\\text\{\\times\}\{10\}^\{\-3\}SCD–2DFM\-PINN \+ first order58,2441\.204×10−31\.204\\text\{\\times\}\{10\}^\{\-3\}5\.916×10−35\.916\\text\{\\times\}\{10\}^\{\-3\}SCD–2DSPINN42,3681\.726×10−31\.726\\text\{\\times\}\{10\}^\{\-3\}4\.543×10−34\.543\\text\{\\times\}\{10\}^\{\-3\}SCD–4DFM\-PINN58,4971\.458×10−21\.458\\text\{\\times\}\{10\}^\{\-2\}1\.440×10−21\.440\\text\{\\times\}\{10\}^\{\-2\}SCD–4DFM\-PINN \+ first order58,5021\.488×10−21\.488\\text\{\\times\}\{10\}^\{\-2\}1\.695×10−21\.695\\text\{\\times\}\{10\}^\{\-2\}SCD–4DSPINN84,7362\.969×10−32\.969\\text\{\\times\}\{10\}^\{\-3\}3\.379×10−33\.379\\text\{\\times\}\{10\}^\{\-3\}SRD–2DFM\-PINN58,2415\.672×10−45\.672\\text\{\\times\}\{10\}^\{\-4\}2\.445×10−32\.445\\text\{\\times\}\{10\}^\{\-3\}SRD–2DFM\-PINN \+ first order58,2443\.242×10−43\.242\\text\{\\times\}\{10\}^\{\-4\}1\.529×10−31\.529\\text\{\\times\}\{10\}^\{\-3\}SRD–2DSPINN42,3681\.629×10−41\.629\\text\{\\times\}\{10\}^\{\-4\}2\.629×10−42\.629\\text\{\\times\}\{10\}^\{\-4\}SRD–4DFM\-PINN58,4971\.272×10−21\.272\\text\{\\times\}\{10\}^\{\-2\}1\.616×10−21\.616\\text\{\\times\}\{10\}^\{\-2\}SRD–4DFM\-PINN \+ first order58,5021\.266×10−21\.266\\text\{\\times\}\{10\}^\{\-2\}1\.390×10−21\.390\\text\{\\times\}\{10\}^\{\-2\}SRD–4DSPINN84,7361\.096×10−31\.096\\text\{\\times\}\{10\}^\{\-3\}9\.514×10−49\.514\\text\{\\times\}\{10\}^\{\-4\}QCD–2DPINN58,1777\.079×10−47\.079\\text\{\\times\}\{10\}^\{\-4\}4\.185×10−34\.185\\text\{\\times\}\{10\}^\{\-3\}QCD–2DFM\-PINN58,2416\.206×10−46\.206\\text\{\\times\}\{10\}^\{\-4\}1\.863×10−31\.863\\text\{\\times\}\{10\}^\{\-3\}QCD–2DFM\-PINN \+ first order58,2447\.652×10−47\.652\\text\{\\times\}\{10\}^\{\-4\}6\.360×10−46\.360\\text\{\\times\}\{10\}^\{\-4\}QCD–2DSPINN42,3682\.833×10−32\.833\\text\{\\times\}\{10\}^\{\-3\}1\.196×10−21\.196\\text\{\\times\}\{10\}^\{\-2\}QCD–4DPINN58,3051\.586×10−21\.586\\text\{\\times\}\{10\}^\{\-2\}1\.671×10−21\.671\\text\{\\times\}\{10\}^\{\-2\}QCD–4DFM\-PINN58,4979\.593×10−39\.593\\text\{\\times\}\{10\}^\{\-3\}1\.051×10−21\.051\\text\{\\times\}\{10\}^\{\-2\}QCD–4DFM\-PINN \+ first order58,5029\.630×10−39\.630\\text\{\\times\}\{10\}^\{\-3\}9\.740×10−39\.740\\text\{\\times\}\{10\}^\{\-3\}QCD–4DSPINN84,7364\.978×10−34\.978\\text\{\\times\}\{10\}^\{\-3\}5\.098×10−35\.098\\text\{\\times\}\{10\}^\{\-3\}QRD–2DPINN58,1779\.459×10−49\.459\\text\{\\times\}\{10\}^\{\-4\}1\.503×10−31\.503\\text\{\\times\}\{10\}^\{\-3\}QRD–2DFM\-PINN58,2415\.174×10−45\.174\\text\{\\times\}\{10\}^\{\-4\}1\.222×10−31\.222\\text\{\\times\}\{10\}^\{\-3\}QRD–2DFM\-PINN \+ first order58,2444\.948×10−44\.948\\text\{\\times\}\{10\}^\{\-4\}8\.335×10−48\.335\\text\{\\times\}\{10\}^\{\-4\}QRD–2DSPINN42,3686\.818×10−46\.818\\text\{\\times\}\{10\}^\{\-4\}8\.145×10−48\.145\\text\{\\times\}\{10\}^\{\-4\}QRD–4DPINN58,3051\.679×10−21\.679\\text\{\\times\}\{10\}^\{\-2\}2\.138×10−22\.138\\text\{\\times\}\{10\}^\{\-2\}QRD–4DFM\-PINN58,4971\.114×10−21\.114\\text\{\\times\}\{10\}^\{\-2\}1\.225×10−21\.225\\text\{\\times\}\{10\}^\{\-2\}QRD–4DFM\-PINN \+ first order58,5021\.320×10−21\.320\\text\{\\times\}\{10\}^\{\-2\}1\.379×10−21\.379\\text\{\\times\}\{10\}^\{\-2\}QRD–4DSPINN84,7366\.139×10−36\.139\\text\{\\times\}\{10\}^\{\-3\}6\.293×10−36\.293\\text\{\\times\}\{10\}^\{\-3\}
### Analysis and Scope of the Evidence

##### High\-dimensional results\.

SPINN is the most accurate method on both four\-dimensional equations in both solution settings\. This is the principal empirical conclusion of the comparison: for the high\-dimensional smooth problems studied here, the rank\-128 separable architecture remains a strong baseline and achieves the lowest error\. The magnitude of its advantage, however, depends strongly on the solution structure\. Under Setting S, the best\-checkpoint error of SPINN is 4\.91 times lower than FM\-PINN without the explicit first\-order output branch on 4D convection–diffusion and 11\.60 times lower on 4D reaction–diffusion\. Under Setting Q, these factors decrease to 1\.93 and 1\.82, respectively\. Thus, SPINN still performs better in four dimensions, but its margin narrows substantially once non\-axis\-aligned quadratic phases are introduced\.

##### Two\-dimensional sensitivity\.

The two\-dimensional comparison changes more sharply\. In Setting S, FM\-PINN is more accurate on convection–diffusion, whereas SPINN is more accurate on reaction–diffusion\. In Setting Q, FM\-PINN without the explicit first\-order output branch is more accurate on both equations: its best\-checkpoint error is 4\.56 times lower than SPINN for convection–diffusion and 1\.32 times lower for reaction–diffusion\. In particular, the 2D reaction–diffusion ordering reverses between the two exact solutions\. This result shows that conclusions drawn from a single manufactured solution can be fragile, even when the PDE operator, optimizer, sample budget, network scale, and evaluation grid are unchanged\.

##### Interpretation\.

These experiments do not establish that either architecture dominates all nonseparable PDE solutions\. Setting S is unusually favorable to a separable representation because its exact solution is explicitly constructed from at most five separable terms, while SPINN uses rank 128\. Setting Q is not an adversarial destruction of separability, but it does introduce coordinate products inside smooth phase functions, which changes the approximation burden of a fixed\-rank CP representation\. FM\-PINN directly forms cross\-coordinate interaction features and consequently shows a clear advantage in the 2D quadratic\-phase problems\. At the same time, the 4D results demonstrate that rank\-128 SPINN can still approximate these coupled solutions effectively\.

Accordingly, we view the present study as evidence of a solution\-structure\-dependent regime rather than a universal ranking\. A more complete investigation should vary interaction strength and frequency, measure the numerical CP rank of the target solution, sweep the SPINN rank, and aggregate multiple random seeds\. Such a controlled study is left for future work\.