Operator Boosting Produces Pareto-Efficient PDE Surrogates

Abstract: 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.

• The paper introduces Operator Boosting, which sequentially trains tiny FNO, DeepONet, or CNO correction operators on residual fields with validation-selected shrinkage.
• Across 30 dataset–architecture pairs, boosted stacks improve mean relative error in 21 cases, reduce parameters by 72–95%, and Pareto-dominate full baselines on 7 of 10 PDE benchmarks.
• The method is most effective when small operators can represent systematic residual structure, but architecture–PDE mismatches can degrade accuracy, highlighting the need for residual-aware model selection.

Overview and motivation

This paper introduces Operator Boosting, a stagewise residual-learning framework for constructing compact neural-operator surrogates for partial differential equations (PDEs). The motivating observation is that full-size monolithic neural operators—Fourier neural operators (FNOs), DeepONets, convolutional neural operators (CNOs)—can be costly to store, transfer, and evaluate in many-query scientific workflows such as design optimization, uncertainty quantification, control, and edge deployment. Rather than training a large model and compressing it afterward via pruning, quantization, tensorization, or distillation, the paper asks whether compact surrogates can be built directly as additive stacks of tiny same-family operators trained sequentially on residual fields.

The framework is a functional-gradient analogue of classical boosting: under squared loss in normalized output coordinates, the negative functional gradient at each stage is the residual field, so each stage reduces to fitting a tiny operator to the current ensemble's error. This connects the method to Friedman's gradient boosting machine, $L_2$-boosting, and early-stopping analyses of boosting, while transferring the residual principle from scalar regression to operator-valued learning between function spaces. A related spectral analysis of FNOs had previously suggested that residual correction can recover information missed by an initial operator; this work systematizes that idea into a model-family-agnostic wrapper evaluated at scale.

Method

The surrogate is represented as an additive operator-valued model

$G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$

where each $H_m^f$ is a tiny neural operator from a fixed family $$f \in \{\mathrm{FNO}, \mathrm{DeepONet}, \mathrm{CNO}\}$$ and $\eta_m \geq 0$ is a shrinkage coefficient selected on validation data from a finite grid $\Lambda$ that includes $\eta = 0$. The initializer $G_0$ is the zero predictor in normalized coordinates (the empirical mean field after denormalization) and contributes no trainable parameters, so all capacity resides in the correction stages.

Each stage proceeds as follows. Given the current stack $G_{m-1}$, residual targets $r_i^{(m)} = \widetilde{u}_i - G_{m-1}(a_i)$ are computed on the training split. A tiny operator is trained on these residuals with mean-squared error plus an optional Fourier-weighted penalty emphasizing high-frequency residual structure. The shrinkage coefficient is then chosen by minimizing validation error over $G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$0, which makes each stage validation-safe: if a trained correction does not reduce validation error, it can be rejected outright ($G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$1). All reported results use the final $G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$2-stage predictor rather than a post hoc best-stage selection, avoiding test-set leakage into stage choice.

Two design choices deserve emphasis. First, comparisons are strictly within-family: a full FNO is compared against boosted tiny FNOs, a full DeepONet against boosted tiny DeepONets, and so on. This isolates the effect of stagewise residual learning from the effect of changing architecture class. Second, all models are trained under a fixed-budget terminal-iterate protocol, with each boosted stage receiving the same epoch budget as the corresponding full baseline. The authors note explicitly that Operator Boosting reduces deployed parameter count but not necessarily training cost, since $G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$3 residual operators must be trained rather than one monolithic model—a distinction that matters for offline surrogate construction but not for repeated deployment.

Evaluation uses test relative $G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$4 error after denormalization, aggregated over ten random seeds with 95% confidence intervals. Pareto dominance is defined empirically: a boosted stack dominates when it achieves no larger mean relative error with strictly fewer trainable parameters.

Experimental setup

The evaluation spans one-, two-, and three-dimensional benchmarks drawn from PDEBench, APEBench, and The Well: 1D advection, Burgers, and reaction–diffusion; 2D Darcy flow, reaction–diffusion, incompressible Navier–Stokes, shallow-water dynamics, and active matter; and 3D compressible Navier–Stokes and magnetohydrodynamics (MHD). Trajectory datasets are converted to supervised transition pairs using states separated by a time offset, with the offset provided as input. Normalization statistics come exclusively from training splits. Each dataset–architecture pair is run with 10 seeds, three boosting stages per stack, and fixed tiny architectures throughout.

Results

Across 30 completed dataset–architecture pairs, boosted stacks improve mean relative error in 21 of 30 pairs, with 17 pairs having entirely positive 95% confidence intervals, while every boosted stack reduces trainable parameter count by approximately 72–95%. In best-model comparisons, boosted stacks empirically Pareto-dominate the best full-size baseline on 7 of 10 PDE benchmarks.

The strongest results are substantial rather than marginal:

Benchmark	Family	Error reduction	Size reduction
2D Navier–Stokes	FNO	74.8%	94.5%
2D shallow water	CNO	51.9%	74.3%
2D shallow water	DeepONet	45.2%	88.1%
1D advection	FNO	44.2%	86.2%
3D compressible Navier–Stokes	DeepONet	27.6%	72.2%
3D MHD	FNO	4.4%	95.2%

The 2D Navier–Stokes result is the headline case: the boosted FNO stack reduces mean relative error from 0.0805 to 0.0188 with a confidence interval of [69.2, 80.4] and 10/10 seed wins, using only about 5.5% of the baseline's parameters. Three-dimensional results support generality beyond low-dimensional problems, though the 3D compressible Navier–Stokes gain for FNO has a confidence interval spanning zero ([−1.8, 36.5]), so it should be read as suggestive rather than conclusive at the within-family level.

The failures are equally informative and are not isolated noise. CNO on 2D Navier–Stokes degrades severely (mean performance delta of −202%, 0/10 wins), FNO on 2D reaction–diffusion loses 20.6% on average, and CNO shows consistent negative deltas on Darcy flow and MHD. Active matter is a failure case for both FNO and CNO. These cases indicate that residual boosting fails when the tiny correction class cannot represent the relevant residual structure, when the full baseline already captures the dominant dynamics, or when residual fitting amplifies architecture-specific weaknesses—consistent with prior diagnostic evidence that neural-operator error modes are strongly PDE- and architecture-dependent.

The best-model Pareto analysis permits architecture selection among boosted stacks. Boosted surrogates dominate on advection, Burgers, 1D reaction–diffusion, 2D Navier–Stokes, Darcy, shallow water, and 3D compressible Navier–Stokes; the three tradeoff cases are 2D reaction–diffusion, active matter, and MHD, where compression is achieved but the best boosted stack does not match the best full-size error. Notably, the best boosted shallow-water surrogate comes from the CNO family even though the best full-size baseline is an FNO, showing that cross-family selection can further improve the frontier.

Discussion

The pattern of successes and failures supports a specific interpretation: boosting helps most when the full model's error contains a systematic component representable by the tiny operator class. In those regimes, stagewise residual fitting changes the optimization problem rather than merely compressing it, since successive small operators concentrate capacity on the remaining error. Where the residual geometry is mismatched to the tiny architecture—as appears to be the case for CNO on Navier–Stokes and FNO on 2D reaction–diffusion—the procedure can actively hurt.

For deployment, the practical claim is that trainable parameter count can often be reduced by roughly an order of magnitude without sacrificing accuracy, and frequently while improving it. The author is careful to scope this claim: parameter count proxies storage, model transfer, and deployability, but is not identical to latency or peak memory. Because all correction operators receive the same input field, their forward passes are independent and parallelizable in principle, but no wall-clock latency measurements are reported.

Limitations and open questions

Several limitations are conceded directly in the paper. The evaluation reports only trainable parameter count and relative $G_M(a) = G_0(a) + \sum_{m=1}^{M} \eta_m H_m^f(a),$5 error; model size in bytes, activation memory, latency, throughput, and energy consumption are not measured, despite boosted stacks containing multiple forward modules. The experimental configuration fixes three stages, fixed tiny architectures, and a fixed shrinkage grid without per-PDE tuning, so the results bound what a simple homogeneous recipe achieves rather than what optimized boosting could achieve. Missing baselines include parameter-matched tiny monolithic models, independently trained tiny ensembles, pruning, quantization, distillation, and tensorization—so the Pareto improvements are established against full-size monoliths only, not against the broader efficiency literature. Finally, the mechanism remains architecture- and PDE-dependent, and the paper leaves open how to diagnose whether a remaining residual is spectral, localized, low-rank, multiscale, or boundary-dominated before choosing a correction family. Heterogeneous boosting, where stages draw from different operator families, is proposed as a natural extension but introduces additional stage-selection and training-cost questions that this work deliberately does not address.

Conclusion

Operator Boosting formulates PDE surrogate construction as stagewise functional-gradient descent in operator space, using tiny same-family neural operators fitted to successive residual fields with validation-selected shrinkage. Across FNO, DeepONet, and CNO backbones and benchmarks spanning 1D to 3D, the approach yields positive mean accuracy gains in 21 of 30 dataset–architecture pairs, uniform parameter reductions of 72–95%, and empirical Pareto dominance over the best full-size baseline on 7 of 10 benchmarks—including 2D Navier–Stokes, shallow-water dynamics, and 3D compressible Navier–Stokes. The failures are systematic rather than incidental and delineate the regime where the tiny correction class cannot represent the residual structure left by earlier stages. The study establishes homogeneous residual stacking as a viable direct route to parameter-efficient neural PDE surrogates, contingent on matching the correction class to the residual geometry of the target problem.