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 ∂-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, ∂-NO consistently improves prediction accuracy, with learned orders that reflect a representation of known physics modulated by noise and finite-sample constraints.
Figures & tables
Figure 1: ∂ -NO architecture. 2 2 2 Code: https://github.com/FaresBMehouachi/delNO The input function u is augmented with learnable fractional derivatives Dβku , where each order βk and scale sk is optimized jointly with backbone weights during training. This simple augmentation wraps any neural operator backbone without modification.
Burgers ( ν=0.1 )
Burgers ( ν=0.001 )
KdV
Darcy (Elliptic)
NS2D ( ν=10−5 )
σ=0
σ=0.05
σ=0
σ=0.05
σ=0
σ=0.05
σ=0
σ=0.05
σ=0
σ=0.05
FNO
0.065
0.250
1.138
1.632
1.598
2.255
5.643
5.590
3.479
3.575
∂ -FNO
0.055
0.190
0.837
1.366
1.549
2.195
5.426
5.415
3.428
3.511
Improv.
15.1%
24.1%
26.5%
16.3%
3.0%
2.7%
3.8%
3.1%
1.5%
1.8%
TFNO
0.071
0.209
0.702
1.032
2.312
2.745
5.223
5.206
3.428
3.534
∂ -TFNO
0.062
0.178
0.676
0.917
2.227
2.660
4.843
4.906
3.228
3.368
Table 1: Relative L2 error (%), averaged over 5 seeds, across 5 PDE benchmarks × 4 backbones × {clean, noisy}. Bold = best per (backbone, dataset, noise) cell. ∂ -augmentation improves every cell tested.
Figure 2: Burgers (smooth, ν=0.1 ): Theory vs. experiment. Top: relative L2 error; bottom: learned β∗ . Theorem 5.6 ’s two unconditional predictions are supported: (left) σ↑⇒β∗↓ ; (center) n↑⇒β∗↑ , with β∗ approaching but not reaching m=2 , consistent with the unreachable recovery threshold (Remark 5.7 ). At smoothness τ≈2.2 (comfortably above the Assumption 4.2 threshold s>m+d/2 ), the capacity formula predicts k↑⇒β∗↓ (right), which holds monotonically.
Figure 3: Burgers (near-shock, ν=0.001 ): Theory vs. experiment. Same panels as Figure 2 . The noise prediction ( σ↑⇒β∗↓ ) holds clearly; the sample-size effect ( n↑⇒β∗↑ ) attenuates to a noisy plateau near the smoothness boundary ( τ≈0.8 ), as expected when Assumption 4.2 is marginal. At rough-data smoothness τ≈0.8 (near the Assumption 4.2 smoothness boundary), the capacity formula admits a weakened k -effect; the empirical k -sweep is correspondingly noisier and non-monotone rather than the clean monotone of the smooth case, consistent with this regime change.
Poisson
NS Shear
Wave
Allen–Cahn
Transp. Cont.
Darcy
Airfoil
σ
0
0.05
0
0.05
0
0.05
0
0.05
0
0.05
0
0.05
0
0.05
L2
CNO
1.788
3.664
9.686
12.218
0.590
1.990
2.618
1.980
0.416
0.413
0.763
0.724
1.184
1.130
∂ -CNO
1.053
1.293
7.104
7.746
0.520
0.900
1.632
1.575
0.353
0.381
0.752
0.718
1.137
1.111
Improv.
41.2%
64.7%
26.7%
36.6%
11.8%
54.8%
37.7%
20.4%
15.2%
7.7%
1.4%
0.8%
4.0%
1.7%
L1
CNO
1.200
2.504
8.120
10.080
0.440
1.550
0.870
0.940
0.358
0.380
0.514
0.482
0.446
0.473
∂ -CNO
0.660
0.884
5.910
6.110
0.439
0.692
0.646
0.857
0.341
0.352
0.501
0.457
0.443
0.410
Table 2: ∂ -CNO vs. CNO on seven benchmarks from the suite of Raonic et al. [30] , evaluated under both relative L2 error and median relative L1 error (averaged over 5 seeds). Bold = best per (dataset, noise, metric) cell. ∂ -CNO improves CNO in 14/14 settings under both metrics.
Appendix figures & tables19 assets
Supplementary material from the paper’s appendix.
Appendix
Notation
Meaning
Spaces and functions
Ω⊂Rd
Spatial domain
d
Spatial dimension
Hs(Ω)
Sobolev space of regularity s
Bs(R)
Sobolev ball {u∈Hs:∥u∥Hs≤R}
u,u0
Input function / initial condition
Appendix
Table 3: Summary of notation used throughout the paper.
Backbone
Width
Modes
Layers
Heads
Slices
FNO
64–128
12–16
4
–
–
Transolver
128–256
–
8
8
32–64
Appendix
Table 4: Architecture hyperparameters by backbone.
Dataset
Resolution
Batch
LR
Scheduler
Tin/Tout
Burgers
1024
8
10−3
OneCycle
–
KdV
1024
8
10−3
OneCycle
–
Darcy
85×85
4
10−3
OneCycle
–
Navier-Stokes
64×64
2–20
5×10−4
Step
10/10
Appendix
Table 5: Training hyperparameters by dataset.
Model
Dataset
Baseline
Fixed {1,2}
Learned {β1,β2}
Learned Values
FNO
Burgers
0.098
0.078
0.079
1.36, 2.00
FNO
Navier-Stokes
9.64
9.46
9.18
1.05, 2.00
Transolver
Burgers
0.65
0.63
0.63
0.74, 1.83
Transolver
Darcy
1.10
0.98
0.95
0.67, 1.78
Appendix
Table 6: Ablation: learnable vs fixed derivative orders. Columns 3–5 show relative L2 error (%); column 6 shows the learned (β1,β2) values. Learning β improves over fixed orders on NS and Darcy. The FNO baseline here uses the original Li et al. [1] implementation (rather than the THUML default of Table 1 ), which accounts for the differing baseline values on NS. The Darcy rows further differ from Table 1 in the unit-Gaussian normalization of the data; absolute values are therefore not comparable across the two tables.
Dataset
d
m
Required s
Fitted α
semp
Satisfied?
Burgers ν=0.1
1
2
>2.5
∼7.2
∼3.1
✓
Burgers ν=0.001
1
2
>2.5
∼5.8
∼2.4
marginal
Darcy
2
2
>3
∼9.5
∼3.75
✓
Navier-Stokes
2
2
>3
∼8.2
∼3.1
✓
Appendix
Table 7: Verification of spectral decay assumption across benchmarks.
Model
Burgers ( ν=0.1 )
Burgers ( ν=0.001 )
Darcy
NS ( ν=10−5 )
FNO
95
78
313
63
∂ -FNO
149
124
664
135
Transolver
412
382
1745
2241
∂ -Transolver
517
480
2201
3158
Appendix
Table 8: Training time per 100 epochs (seconds).
Data-only
Physics-informed ( λpde>0 )
Dataset
σ
FNO
∂ -FNO
Δ∂
PINO
∂ -PINO
Δ∂
Burgers ( ν=10−1 )
0
0.07
0.07
7.1%
1.32
0.69
47.8%
Burgers ( ν=10−1 )
0.05
0.41
0.38
9.3%
1.44
0.67
53.2%
Burgers ( ν=10−3 )
0
0.61
0.56
6.9%
0.62
0.57
8.9%
Burgers ( ν=10−3 )
0.05
1.09
1.02
6.4%
1.05
1.01
3.6%
KdV
0
2.50
2.49
0.1%
6.69
6.19
7.5%
Appendix
Table 9: Four-way physics-informed comparison. Relative L2 error (%). Δ∂ is the relative gain from adding ∂ -augmentation to each baseline. ∂ -PINO improves PINO in 8/8 settings.
σ
FNO
∂ -FNO
Random-FNO
Wider-FNO
0
0.000597
0.000560
0.000638
0.000579
0.01
0.000824
0.000733
0.000843
0.000747
0.02
0.001344
0.001110
0.001299
0.001178
0.05
0.004210
0.003357
0.004037
0.003847
0.10
0.008356
0.006704
0.008811
0.009224
0.25
0.018227
0.013698
0.018047
0.018264
Appendix
Table 10: Matched-parameter ablation on Burgers ( ν=0.1 ). Mean relative L2 error across 5 seeds; all four models at parameter count ≈292 K. Bold = best per row. ∂ -FNO wins 32/35 seeds vs FNO ( 91.4% ), 34/35 vs Random-FNO ( 97.1% ), and 27/35 vs Wider-FNO ( 77.1% ).
Gaussian
Uniform
FNO
LocalNO
FNO
LocalNO
σ
BL
∂
Δ∂
BL
∂
Δ∂
BL
∂
Δ∂
BL
∂
Δ∂
Burg. ν=0.1
0.05
0.00259
0.00210
19.1%
0.00261
0.00215
17.7%
0.00138
0.00111
19.9%
0.00137
0.00113
17.4%
0.10
0.00621
0.00459
26.0%
0.00726
0.00518
28.7%
0.00281
0.00206
26.9%
0.00568
0.00348
38.8%
0.15
0.00944
0.00711
24.7%
0.01264
0.00995
21.3%
0.00530
0.00392
26.0%
0.00908
0.00539
40.6%
0.20
0.01270
0.01050
17.4%
0.01925
0.01584
17.7%
0.00770
0.00598
22.3%
0.01075
0.00757
29.6%
Appendix
Table 11: Stress-test: relative L2 error under Gaussian and uniform additive noise ( σ from 0.05 to 0.30 ), two backbones (FNO, LocalNO). Bold = best per (row, noise type, backbone); Δ∂ in green . ∂ -FNO improves the baseline in 60/60 settings; ∂ -LocalNO improves in 60/60 settings; combined 120/120 wins. Peak improvement: 48.5% (Gaussian) and 48.1% (uniform) on smooth Burgers. n=1000 training samples. Δ∂=(BL−∂-Model)/BL×100 .
Figure 4: Resolution sensitivity: FNO vs. ∂ -FNO on Burgers ( ν=0.1 ) across grid resolutions 64 to 2048 . Left: clean inputs. Right: noisy inputs ( σ=0.02 ). Shaded area indicates the ∂ -FNO improvement. Clean improvements are resolution-stable (6–8%); noisy improvements grow at coarser resolutions, from 14% at 2048 up to 39% at 64 : input-side derivative features are most valuable when the network’s Nyquist headroom is tightest.
Figure 5: Learned β clusters across resolutions (Burgers, K=8 learnable channels). Two clusters consistently emerge near m=2 (diffusion, top) and m=1 (convection, bottom). The diffusion cluster is not resolution-invariant: it drifts below 2 at coarser resolutions, an empirical effect of finite stencil fidelity at large grid spacing rather than a prediction of Theorem 5.6 , whose ξmax -dependence would, with all else fixed, imply the opposite drift. The sweep varies the data spectrum and stencil fidelity alongside ξmax and does not isolate this effect. Downsampling, the most common strategy when GPU memory is limited, silently makes exact physics enforcement increasingly suboptimal. ∂ -NO compensates automatically by learning lower orders.
Figure 6: Directional isotropy of learned orders on two 2D isotropic PDEs. Left: Darcy. Right: Navier-Stokes ( ν=10−5 ). Across the K=4 learnable channels, βx (blue) and βy (orange) track each other closely, with channel means βˉx=1.04,βˉy=1.07 on Darcy ( Δ=0.032 ) and βˉx=1.02,βˉy=1.02 on NS ( Δ=0.002 ). The four channels self-organize into a hierarchy spanning [0.15,2.0] , each capturing a different spectral scale. Isotropy is recovered without being prescribed. (NS channel 4 sits at the cap mmax=2 .)
Figure 7: Learned fractional orders on 3D compressible Navier-Stokes, for three backbones (FNO3D, TFNO, Transolver) under clean ( σ=0 , left) and noisy ( σ=0.05 , right) conditions. Each axis ( x,y,z ) retains a channel with order near 2 (Laplacian, the dominant differential operator), while cross-derivative ( βxy,βxz,βyz ) and triple-derivative ( βxyz ) channels collapse to near-identity ( β≈0 ). The network selectively amplifies the physically dominant differential structure without it being prescribed.
FNO3D
TFNO
Transolver
σ=0
σ=0.05
σ=0
σ=0.05
σ=0
σ=0.05
Model
0.2101
0.2168
0.2079
0.2225
0.1990
0.2039
∂ -Model
0.1964
0.2115
0.1971
0.2033
0.1699
0.1695
Improv.
6.5%
2.4%
5.2%
8.6%
14.6%
16.9%
Appendix
Table 12: 3D Compressible Navier-Stokes: test relative L2 error. ∂ -augmentation improves every (backbone, noise) cell; gains widen under noise on Transolver.
Improv. over FNO
Dataset
σ
FNO
∂ -FNO
Conv-FNO
∂ -FNO
Conv-FNO
Burgers ( ν=0.1 )
0
0.0832
0.0738
0.0791
11.3%
4.9%
Burgers ( ν=0.1 )
0.02
0.1161
0.0899
0.1052
22.6%
9.4%
Darcy
0
5.736
5.474
5.687
4.6%
0.9%
Darcy
0.02
5.721
5.487
5.678
4.1%
0.8%
NS ( ν=10−5 )
0
3.506
3.456
4.219
1.5%
−20.3%
Appendix
Table 13: Conv-FNO vs. ∂ -FNO on three benchmarks. Relative L2 error (%); Δ is the relative improvement over the FNO baseline. ∂ -FNO matches or outperforms Conv-FNO with 34–187 × fewer preprocessing parameters (see note below).
Figure 8: Ablation on the number of derivative channels K∈{1,2,4,8,16,32} . Increasing K consistently improves performance with diminishing returns past K=8 . The natural tradeoff is memory and compute cost. Under noise the curve flattens more slowly, indicating that additional channels help more when the regularization burden is larger.
Figure 9: Sensitivity to GL stencil length KGL∈{2,4,8,16,32,64} . Performance stabilizes around KGL=16 under both clean and noisy conditions. At KGL=2 , accuracy degrades because the stencil is too short to resolve non-integer orders. We use KGL=16 throughout the main paper as a balance between accuracy and compute.
Figure 10: Trajectory of learned β values during training (5 noise levels). All K=8 orders stabilize within 200–400 epochs and converge to values consistent with Theorem 5.6 : higher β∗ under clean data, systematically lower β∗ as σ grows. Initial values are β1=1.0 and β2=2.0 ; the optimizer redistributes channels across the spectral hierarchy as training proceeds.
Figure 11: Interaction between noise level σ and number of derivative channels K . Cell values are the relative improvement of ∂ -FNO over FNO (% reduction in relative L2 ) on Burgers ( ν=0.1 ). Gains broadly increase with both σ and K and concentrate in the high-noise, multi-channel corner, though the pattern is not monotone in either variable at these seed counts. The interaction is super-additive: the largest single gain ( +36.4% ) sits at (σ=0.1,K=8) , in the high-noise regime near the diminishing-returns plateau of Figure 8 , supporting the design choice of multiple channels precisely when data is noisy.
Neural operators are widely used as surrogate solution maps for partial differential equations (PDEs), but full-size models can be costly to store, deploy, and evaluate in many-query scientific workflows. This work introduces Operator Boosting, a stagewise residual-learning framework for constructing compact neural-operator surrogates directly, rather than training a large model and compressing it afterward. Starting from the empirical mean predictor in normalized output coordinates, the method trains a sequence of tiny same-family neural operators on residual fields and incorporates each correction through validation-selected shrinkage. We instantiate the framework with Fourier neural operators (FNOs), DeepONets, and convolutional neural operators (CNOs), and compare boosted tiny stacks against full-size monolithic baselines across one-, two-, and three-dimensional PDE benchmarks from PDEBench, APEBench, and The Well. Across 30 dataset-architecture pairs, 21 show positive mean accuracy gains and 17 have positive confidence intervals, while all boosted stacks reduce trainable parameter count by approximately 72-95%. Best-model comparisons show empirical Pareto improvements on 7 of 10 completed PDE benchmarks, including two-dimensional Navier-Stokes, shallow-water dynamics, Darcy flow, one-dimensional transport and reaction systems, and three-dimensional compressible Navier-Stokes. These results show that Operator Boosting often improves the empirical accuracy-parameter Pareto frontier of neural PDE surrogates, while also exposing PDE- and architecture-dependent regimes where residual boosting fails to offset compression.
Lennon J. Shikhman
Georgia Institute of Technology College of Computing Atlanta, Georgia, USA
Neural operators have emerged as promising surrogate models for solving partial differential equations (PDEs), but struggle to generalise beyond training distributions and are often constrained to a fixed temporal discretisation. This work introduces a physics-informed training framework that addresses these limitations by decomposing PDEs using operator splitting methods, training separate neural operators to learn individual non-linear physical operators while approximating linear operators with fixed finite-difference convolutions. This modular mixture-of-experts architecture enables generalisation to novel physical regimes by explicitly encoding the underlying operator structure. We formulate the modelling task as a neural ordinary differential equation (ODE) where these learned operators constitute the right-hand side, enabling continuous-in-time predictions through standard ODE solvers and implicitly enforcing PDE constraints. Demonstrated on incompressible and compressible Navier--Stokes equations, our approach achieves better convergence and superior performance when generalising to unseen physics. The method remains parameter-efficient, enabling temporal extrapolation beyond training horizons, and provides interpretable components whose behaviour can be verified against known physics.
Vignesh Gopakumar, Ander Gray, Dan Giles +5
UCL Centre for Artificial Intelligence · UK Atomic Energy Authority · LIX, CNRS, ´Ecole Polytechnique +2
Partial differential equations (PDEs) govern complex systems, yet neural operators often struggle to efficiently capture the long-range, nonlocal interactions inherent in their solution maps. We introduce Spectral Filtering Operator (SFO), a neural operator that parameterizes integral kernels using the Universal Spectral Basis (USB), a fixed, global orthonormal basis derived from the eigenmodes of the Hilbert matrix in spectral filtering theory. Motivated by our theoretical finding that the discrete Green's functions of shift-invariant PDE discretizations exhibit spatial Linear Dynamical System (LDS) structure, we prove that these kernels admit compact approximations in the USB. By learning only the spectral coefficients of rapidly decaying eigenvalues, SFO achieves a highly efficient representation. Across six benchmarks, including reaction-diffusion, fluid dynamics, and 3D electromagnetics, SFO achieves state-of-the-art accuracy, reducing error by up to 40% relative to strong baselines while using substantially fewer parameters.
Noam Koren, Rafael Moschopoulos, Kira Radinsky +1
Department of Computer Science Technion - Israel Institute of Technology · Department of Computer Science Princeton University