Fractional is Better: Learnable Derivative Orders in Neural Operator Learning

AI summary

Overview

• Research area: Scientific machine learning — specifically neural operators for solving partial differential equations (PDEs), with a focus on fractional calculus and spectral learning theory.
• Technical level: Advanced. The paper combines Picard iteration on mild solutions, Sobolev-space approximation theory, and a closed-form spectral bias-variance analysis.
• Scope: The paper argues that neural operators should be fed derivative features alongside raw function values, proves that the best derivative order is a learnable fractional order strictly below the PDE order, and validates this on five PDE benchmarks across four backbones plus a convolutional-operator study on seven benchmarks.

What This Paper Is About

Nearly all neural operators take the same input: raw coordinate-value pairs such as (x, u(x)) in 1D or (x, y, u(x, y)) in 2D, ignoring the differential structure that defines the underlying physics. The authors ask whether inputs aligned with PDE differential structure help, and — if so — which derivative orders should be supplied. Their answer is that derivative features do help, but that the optimal derivative order is strictly less than the PDE order at any practically realizable sample size, so the order should be learned from data rather than prescribed.

Key Contributions

1. A derivation of why derivative features help. Using Picard iteration on mild solutions, the authors show that derivatives of the input function appear explicitly in the PDE solution operator, so supplying them reduces the network's implicit differentiation burden.
2. A sharper approximation bound. They prove that the minimax approximation error improves from O(W^{-(s-m)/d}) to O(W^{-s/d}), where W is parameter count, s is input regularity, m is PDE order, and d is spatial dimension. For a second-order PDE in one dimension with s = 3, the exponent improves from 1 to 3.
3. The central finding: β* < m. A closed-form spectral bias-variance analysis shows the MSE-optimal derivative order β* equals m − δ*, with δ* strictly positive whenever n/σ² < e^{A_c}. The gap is an automatic spectral regularizer.
4. ∂-NO (del-NO). An architecture-agnostic augmentation that supplies learnable fractional derivative features to any neural operator backbone, with derivative orders and scales optimized jointly with network weights.

Main Findings

• Optimal order falls short of the physics. Theorem 5.6 gives β* = m − δ* with δ* = (A_c − ln(n/σ²)) / (2 ln(ξ_c · ξ_max)), where A_c := 2s ln ξ_c + (2m + d) ln ξ_max, ξ_c is the network's frequency cutoff, and ξ_max = π/h is the Nyquist frequency. Bias falls and variance rises with β, so the two balance below m.
• Exact physics is astronomically expensive. Recovery of the exact order requires n* = σ² exp(A_c) samples. For a representative FNO setting (unit-interval Burgers, N = 1024, k_max = 16, ξ_max ≈ 3200, ξ_c ≈ 100, A_c ≈ 68), even at σ = 10⁻³ this implies n* ≈ 10²³ samples.
• Three testable monotonicities. The formula predicts σ↑ ⇒ β*↓, n↑ ⇒ β*↑, and k↑ ⇒ β*↓. All three are confirmed empirically; on smooth Burgers (τ ≈ 2.2) the sample-size effect is monotone and β* approaches but never reaches m = 2, while on near-shock Burgers (τ ≈ 0.8) the sample-size effect plateaus and the capacity sweep is noisier and non-monotone, consistent with the smoothness boundary in Assumption 4.2 (s > m + d/2).
• Improvements across every tested cell. Table 1 reports relative L² error (%) averaged over 5 seeds across 5 PDE benchmarks × 4 backbones × {clean, noisy}. The augmentation improves every backbone/dataset/noise combination. Highlights: ∂-FNO improves Burgers (ν = 0.001) from 1.138 to 0.837 clean (26.5%) and from 1.632 to 1.366 noisy (16.3%); ∂-Transolver improves the same benchmark from 18.988 to 4.633 clean (75.6%) and from 7.521 to 4.049 noisy (46.2%); ∂-TFNO improves Darcy from 5.223 to 4.843 (7.3%) and from 5.206 to 4.906 (5.8%); ∂-LocalNO improves NS2D (ν = 10⁻⁵) from 1.994 to 1.984 (0.5%).
• The inductive bias is what matters, not extra capacity. On CNO across seven benchmarks from the suite of Raonic et al., ∂-CNO improves CNO in 14/14 settings under both relative L² and median relative L¹ error. Largest relative L² gains include Poisson (1.788 → 1.053 clean, 41.2%; 3.664 → 1.293 noisy, 64.7%), NS Shear (9.686 → 7.104 clean, 26.7%; 12.218 → 7.746 noisy, 36.6%), Wave (1.990 → 0.900 noisy, 54.8%), and Allen–Cahn (2.618 → 1.632 clean, 37.7%). A matched-parameter control (Appendix L.2) isolates the gains to the derivative inductive bias rather than added capacity.
• Robustness holds under a broader noise sweep. A stress test under Gaussian and uniform additive noise with σ ∈ [0.05, 0.3] confirms robustness across every tested setting (Appendix L.3).
• Sub-physical regularization is the operational regime. The exponential dependence in n* means that for any finite n and σ > 0 with n/σ² < e^{A_c}, using the exact PDE order m is never optimal.
• Distinct from prior fractional-calculus work. The authors position ∂-NO as a fourth role for fractional calculus in deep learning — learnable fractional derivative input features applied as preprocessing — alongside fractional PINNs (fractional operators in the loss), fractional gradient methods (fractional gradients during optimization), and fractional architectures (fractional dynamics inside the network). No modification to equation, optimizer, or architecture is required.

Methodology in Plain English

The authors begin with theory. They study a semilinear PDE written as ∂_t u = ℒu + 𝒩[u], where ℒ is a linear differential operator of order m and 𝒩 is a nonlinearity. Writing the mild solution as an integral equation with Green's function G_t, they run Picard iteration for short times and expand the semigroup as G_s * u_0 ≈ u_0 + sℒu_0 + O(s²). Substituting this into the nonlinearity shows that the solution operator depends explicitly on u_0 and its derivatives up to order m — for Burgers' equation, the features {u_0, ∂_x u_0, ∂_xx u_0, u_0 · ∂_x u_0}. The authors focus on the pure derivative terms and leave the nonlinear product to the network, so their method requires only the PDE order m, not equation-specific terms.

They then formalize this as a feature map Φ[u] = (u, D¹u, …, D^m u), prove it is bi-Lipschitz on a Sobolev ball B_s(R) with upper constant √(1+m), and apply Yarotsky's approximation theory to show the minimax rate improves from W^{-(s-m)/d} to W^{-s/d}.

To decide which order is best, they turn to spectral analysis. Assuming a power spectral density S(ξ) = C_S|ξ|^{-2s-d} (the worst-case envelope for H^s integrability, at which Matérn-s fields sit), they decompose MSE into a bias term penalizing the missing high-frequency content |ξ|^{m-β} above the network's cutoff ξ_c, and a variance term capturing noise amplified by |ξ|^{β} when computing D^βu. Both integrals have closed forms; minimizing their sum yields the closed-form β*.

For implementation, they define the Weyl fractional derivative via the Fourier multiplier (iξ)^β = |ξ|^β e^{i(π/2)β sgn(ξ)} with the PDE-agnostic derivative ∂/∂β (iξ)^β = (iξ)^β (log|ξ| + iπ/2 sgn(ξ)), which makes gradient-based learning of β possible. Because Weyl derivatives can produce Gibbs ringing near discontinuities or on non-periodic domains, the implementation uses the Grünwald-Letnikov finite-difference form with generalized binomial coefficients extended to non-integer β via the Gamma function. The h^{-β} prefactor is dropped in practice for training stability, since the learnable scale s_k and the network absorb any constant scaling.

Concretely, ∂-NO replaces the input (x, u) with (x, u, s_1𝔇^{β_1}u, s_2𝔇^{β_2}u, …, s_k𝔇^{β_k}u), where both the orders β and the scales s are learned jointly with the backbone weights. For the second-order benchmarks (Burgers, Darcy, NS) with m = 2 they use two features initialized at β_1 = 1.0 and β_2 = 2.0; for third-order KdV (m = 3) a third feature initialized at β_3 = 3.0. For 2D problems they add directional derivatives along x and y plus one mixed partial with two learnable orders and a single shared scaling factor. Models train with AdamW, with the β parameters on a 10× higher learning rate (10⁻² vs 10⁻³) because they sit at the input level and their gradients attenuate through the full backbone.

Why This Matters

Impact on research. The paper challenges the intuition that surrogate models should be handed exactly the physics they are meant to approximate. It shows that matching the true PDE order is provably suboptimal at finite sample size, and that learning the derivative order from data acts as automatic spectral regularization. It also opens a bridge between fractional calculus and operator learning that is orthogonal to existing work — the augmentation wraps any backbone without touching the equation, the optimizer, or the architecture. Because the analysis is spectral, it yields falsifiable predictions about how the learned order should shift with noise, data volume, and network capacity, and those predictions are tested experimentally rather than asserted.

Real-world applications. The introduction names climate modeling, materials science, and engineering design as domains where neural operators serve as fast PDE surrogates. Any setting where a learned surrogate replaces an expensive numerical solver stands to benefit — for example, repeated forward simulations inside a design or uncertainty-quantification loop, real-time control built on a learned flow model, or inverse problems where the forward map is queried thousands of times.

Industry relevance. The method adds modest overhead and requires only the PDE order, not a fully specified equation, which matters in industrial settings where the governing physics is unknown or only partially specified. Its largest average gains appear on attention-based backbones (up to 75.6% relative error reduction on near-shock Burgers with Transolver) and on convolutional operators (14/14 settings improved), suggesting the augmentation is most valuable exactly where the backbone would otherwise have to learn differentiation from finite, noisy samples.

Future Directions

• Learning nonlinear feature products too. The theory shows u_0 · ∂_x u_0 emerges from the convective nonlinearity, but this work deliberately leaves the nonlinear product to the network. Extending the augmentation to learnable derivatives and their nonlinear combinations is a natural next step.
• Pushing theory validation past the smoothness boundary. The sample-size and capacity predictions weaken or become noisy as τ approaches the Assumption 4.2 boundary (τ ≈ 0.8 on near-shock Burgers versus τ ≈ 2.2 on smooth Burgers). A sharper theory in the marginal-regularity regime would clarify when the framework applies.
• Non-periodic and irregular domains. The Grünwald-Letnikov formulation is introduced to avoid Fourier artifacts on non-periodic domains; its behavior relative to the Weyl formulation, beyond the training-stability rationale given, is left largely to the appendix.
• Composition with physics-informed training. The paper reports complementarity with PINO physics-informed training (Table 9), and reports robustness analyses including emergent directional isotropy on 2D/3D isotropic PDEs (Figs. 6, 7) and resolution-dependent gains (Fig. 4) in the appendix — each of these is a candidate direction for further systematic study.

Target Audience

This paper is aimed at machine learning and computational science researchers working on neural operators, PDE surrogates, and scientific machine learning, particularly those interested in inductive biases for function-space learning. It will also appeal to readers with a background in fractional calculus or approximation theory who want to see those tools applied to operator learning. Practitioners building surrogate models for simulation-heavy workflows can use the ∂-NO recipe directly — the augmentation is architecture-agnostic, and the accompanying code is released at https://github.com/FaresBMehouachi/delNO — but the paper assumes comfort with Sobolev spaces, Fourier analysis, and bias-variance reasoning.

Authors’ abstract

Neural operators learn mappings between function spaces, enabling fast surrogate solutions to partial differential equations. Despite remarkable architectural diversity, these methods often share a common input representation: raw coordinate-value pairs $(x, u(x))$. We ask whether inputs aligned with PDE differential structure can improve learning. Through Picard iteration on mild solutions, we show that derivatives of the input appear explicitly in the solution operator, suggesting that providing derivative features should reduce the network's implicit differentiation burden. We prove this intuition: providing derivative features improves approximation bounds from $O(W^{-(s-m)/d})$ to $O(W^{-s/d})$, where $W$ is parameter count, $s$ is input regularity, $m$ is the PDE order, and $d$ is spatial dimension. Our central finding, however, is a surprise: the optimal derivative order $β^*$ is strictly less than the PDE order $m$. This gap arises from a bias-variance tradeoff in spectral space that we characterize in closed form. Learning $β$ from data achieves automatic spectral regularization. We introduce $\partial$-NO (del-NO), for derivative-augmented neural operators, an architecture-agnostic augmentation that provides learnable fractional derivative features to any neural operator backbone. Across benchmark problems and architectures, $\partial$-NO consistently improves prediction accuracy, with learned orders that reflect a representation of known physics modulated by noise and finite-sample constraints.