Fractional is Better: Learnable Derivative Orders in Neural Operator Learning

• Raw-value inputs omit PDE differential structure, leaving neural operators to learn differentiation implicitly from finite, noisy samples.The paper asks whether inputs aligned with that structure can improve learning.Source: §1 · pp. 1–2
• For practically realizable finite, noisy datasets, the optimal derivative order β* is strictly below the PDE order m.The gap shrinks as data grows, but remains positive at every practically realizable scale.Source: §1 · p. 2
• Because the optimal order depends on unknown noise and sample size, the paper proposes learning it directly from data.The resulting features balance physical fidelity against statistical efficiency.Source: §1 · p. 2
• ∂-NO augments any neural-operator backbone with fractional derivative features whose orders are jointly optimized during training.It wraps existing architectures without modification; the paper distinguishes this preprocessing role from fractional losses, optimizers, and architectures.Source: §1 · p. 2
• ∂-NO consistently improves prediction accuracy across benchmark problems and backbone choices.The learned orders also reflect known physics modulated by noise and finite-sample constraints.Source: §1 · p. 2

2.1 Neural Operators

• A neural operator Gθ is a finite-dimensional approximation of the true solution operator G, learned from paired input-output examples.The target operator typically maps initial or boundary conditions to PDE solutions.Source: §2.1 · p. 3

2.2 A Concrete Example: Burgers’ Equation

• For Burgers’ equation, diffusion and convection make the solution operator depend on second- and first-order derivatives, which standard operators must learn from data.This motivates learning derivative orders rather than prescribing PDE-matched orders, especially when physics is unknown or only partly specified.Source: §2.2 · p. 3

2.3 Fractional Derivatives

• The Weyl fractional derivative supplies a continuous order β: β=1 recovers ∂x, β=2 recovers ∂xx, and intermediate values interpolate smoothly.Its differentiability with respect to β enables gradient-based learning of derivative orders.Source: §2.3 · p. 3

3.1 Mild Solutions and the Green’s Function

• For a semilinear PDE with an order-m linear operator and polynomial or quasilinear nonlinearity, the mild solution is expressed using the Green’s-function semigroup.For constant-coefficient operators, the semigroup has the explicit form Gt=e^tL.Source: §3.1 · p. 4

3.2 The Picard Mechanism

• Picard iteration starts from the semigroup applied to the initial condition; its first iterate exposes derivative structure through the expansion of the nonlinearity.For short times, Gs*u0 is approximated by u0+sLu0+O(s^2), and the nonlinearity expands using its Fréchet derivative.Source: §3.2 · p. 4
• The solution operator is approximated by GT*u0 plus T times an averaged-Green’s-function convolution with N[u0], with a remainder bounded by C(R)T^2.The nonlinearity contains derivatives of u0 up to order m and powers up to its degree.Source: §3.2 · p. 4
• Derivatives enter through both the linear operator’s direct action and nonlinear combinations, and this structure extends beyond short-time parabolic evolution.For Burgers’ equation, the features include u0, ∂xu0, and ∂xxu0, while the nonlinear product u0·∂xu0 is left to the network; the cited discussion also covers stationary, quasilinear, and other evolution equations.Source: §3.2 · p. 4

4 Approximation Theory: Why Features Help

• Derivative features improve the approximation bound by aligning inputs with PDE structure, while selecting the optimal derivative order is deferred to Section 5.The two questions are distinct: whether derivative features help, and which derivative order is best.Source: §4 · p. 5

4.1 The Feature Map

• The feature map pairs u with derivative features; under s > m + d/2, it is bi-Lipschitz on the stated Sobolev ball.The upper stability constant √(1 + m) is independent of the ball radius.Source: §4.1 · p. 5

4.2 Main Approximation Result

• Full derivative features improve the approximation rate from W^(-(s-m)/d) to W^(-s/d) for a PDE solution operator of order m.For a second-order PDE in one dimension with s = 3, the exponent improves from 1 to 3.Source: §4.2 · pp. 5–6

5.1 The Core Tension

• Using β < m leaves high-frequency information uncaptured, while larger β amplifies high-frequency noise; the optimal order balances these opposing effects.The analysis characterizes this as an operator-level bias-variance tradeoff.Source: §5.1 · p. 6, §5 · p. 6

5.2 Spectral Analysis

• Under the stated spectral-decay assumption, bias decreases and variance increases with β, yielding a unique β* that minimizes total MSE.The variance monotonicity condition is ξmax > e^(1/(2β+d)), which the text says holds for typical discretizations with ξmax ≫ 1.Source: §5.2 · pp. 6–7

5.3 The Optimal Derivative Order

• The optimal gap depends on data quality and network capacity: improving signal-to-noise ratio moves β* toward m, while higher frequency capacity can favor lower β*.The capacity effect is less sharp near the smoothness boundary.Source: §5.3 · p. 7
• For finite n and positive noise when n/σ^2 < e^Ac, the optimal order is strictly below m; as data increase or noise vanishes, β* approaches m.The exact-physics recovery threshold is described as practically unreachable at representative settings.Source: §5.3 · p. 7

6 ∂-NO: Derivative-Augmented Neural Operators

• ∂-NO adds fractional derivative features to an existing neural-operator backbone without modifying the backbone architecture.The derivative orders and scales are jointly learned with network weights.Source: §6 · pp. 7–8
• The Grunwald-Letnikov finite-difference formulation supports non-integer orders and is differentiable with respect to β, enabling gradient-based order optimization.It is used instead of the Weyl formulation to avoid Fourier artifacts near discontinuities.Source: §6 · p. 8
• The implementation omits h^-β because learnable scales and the network absorb constant scaling, while the omission helps stabilize training.Augmented backbones are denoted with the ∂-prefix, such as ∂-FNO.Source: §6 · p. 8

7.1 Setup

• The experiments use Burgers, Darcy flow, Navier-Stokes, and KdV datasets with FNO, TFNO, LocalNO, and Transolver backbones.The setup uses learnable derivative features and reports relative L2 errors.Source: §7.1 · p. 8

7.2 Validation of ∂-NO Improvements

• ∂-augmentation improves every tested Table 1 backbone, dataset, and noise combination, with relative L2 errors averaged over five seeds.Table 1 covers five PDE benchmarks, four backbones, and clean and noisy regimes.Source: §7.2 · p. 8, Table 1 · p. 9

7.3 Theory Validation

• On Burgers, learned β* falls as noise rises; with more samples it rises toward, but does not reach, m = 2 on smooth data.The sample-size effect plateaus on near-shock data, and the capacity trend is less clean near the smoothness boundary.Source: §7.3 · pp. 9–10, Figure 2 · p. 10, Figure 3 · p. 10

7.4 Differential Inductive Bias: ∂-CNO

• ∂-CNO improves CNO in 14/14 settings across both relative L2 error and median relative L1 error.Reported peak gains reach 55–65% on Poisson and Wave and approximately 39% on NS Shear.Source: §7.4 · pp. 10–11, Table 2 · p. 11

7.5 Robustness and Directionality

• Learned orders show emergent isotropy on isotropic PDEs, while derivative-feature gains are larger at coarser resolutions and complement PINO in all 8/8 tested settings.On 3D compressible Navier-Stokes, axis-aligned channels settle near β ≈ 2 while the cross-derivative channel approaches identity.Source: §7.5 · p. 11

Neural operators

• ∂-NO wraps neural-operator backbones with learnable derivative features without modifying their architectures.The section positions it as orthogonal to developments including spectral convolutions, attention, and implicit representations.Source: §8 · p. 11

Local-feature operators

• Unlike local-feature operators that add convolutions or localized kernels, ∂-NO supplies fractional differentiation as a local operation with a PDE-structure inductive bias.∂-CNO experiments are described as isolating this bias from generic local processing.Source: §8 · p. 11

Derivative and physics-aware features

• ∂-NO differs from equation-aware and physics-informed methods by learning derivative orders and supplying derivatives as input features rather than prescribing PDE terms or using derivative loss terms.The section also distinguishes its preprocessing role from fractional PINNs, fractional-gradient optimization, and fractional architectures.Source: §8 · pp. 11–12

Spectral perspectives

• The analysis uses neural operators’ spectral structure: fixed Fourier-mode capacity gives a sharply defined cutoff, while spectral bias motivates derivative preprocessing for high-frequency content.Fractional derivatives are treated as learnable feature transformations, not as physical models of anomalous diffusion or long-range dependence.Source: §8 · p. 12

Discussion

• Providing derivative features improves approximation bounds from O(W^-(s-m)/d) to O(W^-s/d), with W parameters, input regularity s, PDE order m, and dimension d.The authors attribute this gain to making derivatives that enter the solution operator directly available as features.Source: §9 · p. 12
• The optimal derivative order β* is less than PDE order m, with the gap characterized by a spectral bias-variance tradeoff and a closed-form expression responsive to noise, sample size, and capacity.For fixed-Fourier-mode operators, the capacity cutoff is sharply defined; learned β is reported to provide automatic spectral regularization.Source: §9 · p. 12
• The sharp predictions of Theorem 5.6 require s > m + d/2 and become directional guidance near shocks; the theorem characterizes population-optimal β*, not the optimizer’s converged value.The method assumes a regular grid for its Grünwald-Letnikov stencil, preprocesses only spatial inputs, and lacks formal trajectory guarantees.Source: §9 · p. 12
• The authors argue that at finite, noisy scales, parameterized rather than prescribed physics lets the optimizer find a compromise between approximation power and estimation variance.They conclude that statistically optimal inductive biases need not match the literal governing equation.Source: §9 · p. 12

Impact Statement

• The impact statement notes potential applications across scientific fields and reports no specific ethical concerns.Source: Impact Statement · p. 13

A Table of Notation

• Table 3 is identified as a summary of notation and its meanings.Source: Table 3 · p. 17

B Proof of Theorem 3.1 (Derivative Emergence)

• The proof establishes that the solution operator depends on derivatives of the initial condition through Picard iteration.For a linear differential operator of order m, the proof identifies dependence on derivatives through order m, along with powers arising from polynomial nonlinearities.Source: B Proof of Theorem 3.1 (Derivative… · pp. 18–19

B.1 Explicit Features for Specific PDEs

• The PDE examples include Burgers, 2D-vorticity Navier–Stokes, Darcy flow, and KdV, with the KdV feature set including derivatives through order three.Source: §B.1 · p. 19

B.2 Extension to Arbitrary Times and Stationary Problems

• The derivative-dependence structure is extended to arbitrary times, stationary problems, and quasilinear systems, where providing input derivatives gives the network direct access to needed quantities.For Navier–Stokes, the bounded inverse Laplacian does not introduce new dependence on initial-condition derivatives.Source: §B.2 · pp. 19–20

B.3 General Evolution Equations: A Spectral Argument

• A spectral argument extends derivative structure to general evolution equations, including hyperbolic and fractional-in-time cases; spatial order m determines the feature structure, not temporal order a.The transfer claim is conditioned on matching the spatial order; the operator norm also depends on a and T.Source: §B.3 · pp. 20–21

C Proof of Proposition 4.3: Bi-Lipschitz Stability

• The stability proof establishes a lower bound with constant 1 because the augmented feature map retains the original function as its first component.It also addresses the upper Lipschitz bound using Sobolev norms.Source: C Proof of Proposition 4.3:… · p. 21

D Proof of Theorem 4.4: Approximation Rates

• Derivative features improve the established approximation bound from O(W^(-(s-m)/d)) to O(W^(-s/d)), but the bounds use different input-side metrics and do not prove a separation.The comparison records the achievable bounds established with and without features; under a common metric, the exponents coincide.Source: D Proof of Theorem 4.4: Approximation… · pp. 22–23

E.1 Assumptions

• The analysis assumes bandlimited inputs and fixed network capacity: derivative preprocessing changes the input representation, not the network’s representational capacity.The capacity cutoff ξc is independent of input features; in FNO-style architectures, retained modes still impose a hard frequency cutoff.Source: §E.1 · p. 23

E.2 Bias Derivation

• The bias–variance tradeoff arises because sub-order features leave the network to learn the spectral mapping from D^βu to D^mu, while higher β amplifies observation noise.The residual requires a frequency multiplier |ξ|^(m-β); noise power in the derivative features scales as σ^2|ξ|^(2β).Source: §E.2 · p. 23, §E.3 · p. 24, §F.2 · p. 26

E.4 Monotonicity and Existence of Optimum

• The total error has a unique global minimum because bias decreases while variance increases with derivative order, producing an interior optimum in the intermediate data regime.The proof establishes existence for 0 < β* < m in that regime and strict convexity of MSE(β) = B(β) + V(β).Source: §E.4 · p. 25

F.2 Derivation

• When n/σ^2 < e^Ac, the optimal order is strictly below m; as n increases or noise vanishes, β* approaches m, and for n/σ^2 ≥ e^Ac no regularization is needed.To leading order, the gap δ* decreases linearly with ln(n/σ^2); the result requires s > m.Source: §F.2 · p. 26, §F.3 · p. 26

F.4 Beyond Leading Order: The O(1/ ln ξc) Correction

• The next-order correction to the optimal gap is O(1/ln ξc) when the interior expansion applies, and is below learned-β measurement resolution at benchmark-scale cutoffs.The expansion applies when 0 < δ* < min(m, s); outside that range the optimum is at an admissible boundary.Source: §F.4 · p. 27

G Spectral Regularization and Tikhonov Interpretation

• Using β < m acts as spectral regularization: it preserves low-frequency signal while increasingly suppressing high-frequency noise, with learned β tuning this effect to the problem’s signal-to-noise ratio.The appendix describes this as a Tikhonov interpretation of sub-physical derivative orders.Source: G Spectral Regularization and Tikhonov… · p. 27
• For β < m, the effective penalty vanishes at low frequencies and grows as |ξ|^{2δ} at high frequencies, where δ = m − β.Relative to using order m directly, the penalty ratio grows as |ξ|^{2(m−β)}.Source: G Spectral Regularization and Tikhonov… · p. 28

H Gr¨unwald-Letnikov Implementation

• For non-periodic domains, the centered Grunwald-Letnikov implementation costs O(N KGL) with KGL = 16 and avoids FFT boundary artifacts.At integer orders, it recovers standard finite differences; its weights vary smoothly with β, enabling gradient-based order optimization.Source: H Gr¨unwald-Letnikov Implementation · p. 28

I Experimental Details

• Experiments generally use 1,000 training and 200 test samples, 500 AdamW epochs, weight decay 10^-5, gradient clipping at 1, and five random seeds.The sample counts can differ for the sample-size sweep.Source: I Experimental Details · p. 28

I.1 Architecture Hyperparameters

• Backbone widths vary by dataset, Transolver uses 64 or 32 slices, and FNO uses 16 Fourier modes for Burgers and 12 for Darcy, NS, and KdV.FNO widths are 64 for Burgers and NS or 128 for Darcy; Transolver widths are 128 for Burgers and Darcy or 256 for NS.Source: §I.1 · p. 29, Table 4 · p. 29

I.2 Training Hyperparameters

• Training uses OneCycle pct start = 0.3, a step scheduler with size 100 and γ = 0.5, and NS batch sizes of 20 for FNO or 2 for Transolver.Dataset-specific training hyperparameters are also summarized in Table 5.Source: §I.2 · p. 29, Table 5 · p. 29
• The KdV dataset maps initial states to time T = 1 using 1,000 training and 200 test trajectories, with pseudo-spectral integrating-factor integration at Δt = 10^-4.Data use δ = 2.5 × 10^-3 on a periodic [0, 2π) domain discretized at N = 1024.Source: §I.2 · p. 29

I.3 ∂-NO Augmentation

• ∂-NO uses GL stencil length 16, maximum derivative order 2 (3 for KdV), and β learning rate 10^-2, with orders initialized at 1 and 2 (plus 3 for KdV).For 2D problems, five derivative features use six learnable orders, including two independently learned orders for the mixed partial.Source: §I.3 · p. 29

I.4 Evaluation

• Evaluation reports test-set-averaged relative L2 error, with confidence intervals computed across five random seeds.Source: §I.4 · p. 29

I.5 Ablation: Learnable vs Fixed Derivative Orders

• Learning β improves over fixed orders on Navier-Stokes and Darcy, but provides no benefit over the fixed-{1, 2} baseline on simpler Burgers benchmarks.The lower learned order β1 falls below the PDE order m = 2, while β2 sits at or near the cap mmax = 2.Source: §I.5 · pp. 29–30, Table 6 · p. 30

J.1 Method

• The verification procedure estimates input smoothness by fitting a power law to empirical power spectra and checks whether it exceeds m + d/2.The fit uses log-log axes over a frequency range excluding DC and Nyquist artifacts.Source: J Verifying the Spectral Decay… · p. 30, §J.1 · p. 30

J.2 Results

• The spectral-decay check finds that most benchmarks satisfy the regularity assumption, but low-viscosity Burgers (ν = 0.001) falls short by approximately 4%.Its estimated smoothness is about 2.4 against the requirement s > 2.5; performance remains good, so the bounds may not be tight and theory is qualitative in this case.Source: §J.2 · pp. 30–31

K Computational Overhead

• Reported training-time overhead is 25–41% for Transolver and 57–114% for FNO, with timings serving as upper bounds because the implementation is unoptimized and runs in eager mode.The derivative computation uses a local 16-term stencil and parallelizes across batch and spatial dimensions; attention-heavy Transolver has lower relative overhead.Source: K Computational Overhead · p. 31, Table 8 · p. 31

L.1 Physics-Informed Comparison: ∂-PINO

• ∂-augmentation improves PINO in all 8/8 settings, while ∂-FNO outperforms PINO across the benchmark despite using no PDE-form information at training time.The two approaches are complementary: adding derivative features also improves PINO.Source: §L.1 · pp. 31–32, Table 9 · p. 31

L.2 Matched-Parameter Ablation

• At matched parameter counts, ∂-FNO wins 32/35 seeds against FNO, 34/35 against Random-FNO, and 27/35 against Wider-FNO.The controls indicate that neither arbitrary local filters nor extra input-layer parameters alone account for the improvement.Source: Table 10 · p. 32, §L.2 · p. 32

L.3 Stress Tests under Gaussian and Uniform Noise

• ∂-augmentation improves all 120 tested settings across five PDE benchmarks, two backbones, two noise distributions, and six noise levels.Peak gains on smooth Burgers at σ=0.30 are 48.5% for Gaussian noise and 48.1% for uniform noise.Source: §L.3 · p. 32, Table 11 · p. 33

L.4 Resolution Sensitivity / KdV

• Clean-input improvements remain 6–8% across resolutions, while noisy-input improvements rise from 14% at resolution 2048 to 39% at 64.At coarser resolutions, ∂-NO learns lower orders, whereas fixed-order approaches cannot compensate in the same way.Source: Figure 4 · p. 33, KdV · p. 33

L.5 Directional Isotropy on 2D Isotropic PDEs

• On Darcy and 2D Navier-Stokes, learned βx and βy agree within Δ < 0.035 without prescribed isotropy, while channels span a hierarchy of spectral scales.The channel hierarchy indicates that the four-channel budget is used across low through high frequencies rather than collapsing to identical orders.Source: §L.5 · p. 34, Figure 6 · p. 34

L.6–L.7 Three-Dimensional Compressible Navier-Stokes

• For 3D compressible Navier-Stokes, each axis retains an order near the Laplacian order m=2, while cross- and triple-derivative channels approach identity across three backbones and noise conditions.Quantitative test errors improve in every backbone/noise cell, with Transolver gains largest both clean and noisy.Source: §L.6 · p. 35, Table 12 · p. 35, §L.7 · p. 35

L.8 Conv-FNO Standalone Comparison

• ∂-FNO matches or outperforms Conv-FNO with 34–187× fewer preprocessing parameters.The comparison covers Burgers, Darcy, and Navier-Stokes.Source: §L.8 · p. 35, Table 13 · p. 36

L.9 Ablations on K and KGL

• Increasing derivative channels improves performance with diminishing returns past K=8, while performance stabilizes around KGL=16 and degrades when KGL=2.More channels cost proportionally more memory; shorter stencils save compute but degrade fractional-order resolution.Source: Figure 8 · p. 36, Figure 9 · p. 36, §L.9 · p. 36

L.10–L.11 Learned-Order Trajectory and Noise–Channel Interaction

• The eight learned orders stabilize within 200–400 epochs, shift lower as noise increases, and yield a largest single joint noise/channel gain of +36.4% at σ=0.1, K=8.The heatmap interaction is superadditive but not monotone in either variable at these seed counts; learned orders redistribute across channels rather than collapsing to redundant values.Source: §L.10 · p. 37, Figure 10 · p. 37, Figure 11 · p. 37