The Symbol of the Surrogate: Measuring Numerical Provenance in Neural PDE Solvers

Ridham Patel

cs.LGOct 7, 2026

The Symbol of the Surrogate: Measuring Numerical Provenance in Neural PDE Solvers

Authors: Ridham Patel

Organizations: Computer Science & Engineering Indian Institute of Technology Gandhinagar Gandhinagar, India

Abstract

Neural PDE surrogates are trained on numerical solver outputs that contain both physical evolution and solver-specific discretization errors. Because surrogates are also evaluated against held-out trajectories from the same solver, standard benchmarks cannot distinguish fidelity to the exact evolution from imitation of the numerical scheme. We introduce an empirical Fourier-symbol diagnostic that probes a trained surrogate's linearized one-step operator with individual Fourier modes and compares it with both exact-evolution and training-scheme references. To address architectural spectral bias, we train identical networks on schemes with orthogonal dissipative and dispersive signatures and compare their learned operators. In linear advection, the learned surrogates reproduce the training schemes' amplitude and phase errors, with the twin-scheme difference reaching more than 99.8% of the analytically predicted full-imitation ceiling. The same behavior occurs for a non-local Fourier neural operator and at the operator level for nonlinear Burgers dynamics. These results show that agreement with solver-generated test data does not by itself establish fidelity to the exact evolution. Fourier-symbol measurements provide a direct diagnostic of numerical provenance.

Figures & tables

model	learned $\|\hat{G}_{\theta}\|$	scheme $\|G_{\mathrm{scheme}}\|$	learned lag	scheme lag	$\alpha$
	reproduced dissipation		reproduced dispersion
upwind	$0.395\pm 0.000$	0.395	$0.003\pm 0.003$	$0$ (exact)	$1.000\pm 0.000$
Crank–Nicolson	$1.000\pm 0.000$	1.000	$0.824\pm 0.000$	0.824	$1.000\pm 0.000$
exact (control)	$0.987\pm 0.006$	–	$0.025\pm 0.003$	–	–

Table 1: Linear advection: the learned symbol matches the training scheme rather than the exact reference. Values are means over the upper half of the resolved spectrum,

\pm

s.d. over 3 seeds. “Learned

|\hat{G}_{\theta}|

” is the modulus of the measured symbol (exact reference:

1

) and “learned lag” is

\Delta_{\mathrm{disp}}

, the phase error relative to the exact reference (exact reference:

0

); each is paired with the value the training scheme itself takes, so that the deviations in the table are read as scheme values reproduced rather than as surrogate errors. The upwind scheme has exactly zero phase error at

\sigma=1/2

(Sec. 2.3 ), and the Crank–Nicolson scheme is exactly unitary.

\alpha

is the provenance coefficient of ( 6 ); it is undefined for the exact control, which has no scheme to attribute to.

Figure 1: Linear advection, per wavenumber $m$ (single seed; Table 1 aggregates over seeds). (a) modulus of the measured symbol against wavenumber: the upwind-trained surrogate (red dots) lies on the upwind scheme curve (red line) across the spectrum, while the Crank–Nicolson-trained surrogate (blue squares) stays at $|\hat{G}_{\theta}|=1$ like its unitary scheme; the exact-trained control (green triangles) scatters about $1$ at high $m$ . (b) phase error $\Delta_{\mathrm{disp}}$ against wavenumber: only the Crank–Nicolson surrogate departs from zero, and it tracks the Crank–Nicolson scheme curve. (c) the two error coordinates plotted against each other, colour-coded by wavenumber: the upwind surrogate moves along the dissipation axis, the Crank–Nicolson surrogate along the dispersion axis, and the exact control stays near the origin (star). This orthogonal separation is the takeaway: each surrogate deviates from the exact reference in the direction its own training scheme deviates, and the tangential arm is what spectral bias cannot produce. (The Crank–Nicolson arm bends leftward at the highest wavenumbers, where the phase error approaches $1$ radian and the first-order reading $\Delta_{\mathrm{diss}}\approx|\hat{G}_{\theta}|-1$ of ( 5 ) no longer applies.) (d) the provenance coefficient $\alpha(k)\approx 1$ for both surrogates across the resolved spectrum, and the two components of the twin-scheme difference.

	corpus	final MSE	mean $\|\hat{G}_{\theta}\|$	mean lag	$\alpha$
CNN (local)	upwind	$1.2\times 10^{-7}$	$0.395\pm 0.000$	$0.003\pm 0.003$	$1.000\pm 0.000$
	Crank–Nicolson	$3.6\times 10^{-8}$	$1.000\pm 0.000$	$0.824\pm 0.000$	$1.000\pm 0.000$
	exact (control)	$3.3\times 10^{-3}$	$0.987\pm 0.006$	$0.025\pm 0.003$	–
FNO (non-local)	upwind	$2.9\times 10^{-8}$	$0.395\pm 0.000$	$0.004\pm 0.001$	$1.000\pm 0.000$
	Crank–Nicolson	$6.5\times 10^{-9}$	$1.000\pm 0.000$	$0.825\pm 0.000$	$1.001\pm 0.000$
	exact (control)	$3.6\times 10^{-8}$	$0.999\pm 0.000$	$0.001\pm 0.000$	–

Table 2: Architecture control on linear advection (mean

\pm

s.d., 3 seeds). Final MSE is against the corpus the model was trained on;

|\hat{G}_{\theta}|

and lag are the measured symbol statistics of Table 1 ;

\alpha

is the provenance coefficient against the training scheme. The non-local FNO reproduces its training scheme as strongly as the local CNN. The exact-trained control (third row of each block) separates the architectures: the CNN cannot fit the non-local exact evolution, whereas the FNO fits it to

3.6\times 10^{-8}

and lands at the exact reference, which is what makes the FNO’s behaviour on scheme corpora a consequence of the training data rather than of its hypothesis class.

surrogate	$D_{\mathrm{scheme}}$ (vs. own solver)	$D_{\mathrm{exact}}$ (vs. spectral ref.)	ratio	scheme artifact
upwind	0.564	3.02	$5.4\times$	3.48
centered	0.127	3.40	$26.8\times$	3.46

Table 3: Forced Burgers, one-step operator fingerprint evaluated on each surrogate’s own training-scheme steady states.

D_{\mathrm{scheme}}

is the spectral distance between the surrogate’s output and its own solver’s next state;

D_{\mathrm{exact}}

is the distance between the surrogate’s output and the dealiased spectral reference operator applied to the same inputs; “scheme artifact” is the corresponding distance between the two reference operators, i.e. the size of the discretization error available to be imitated. The learned operator is

5

–

27\times

closer to its training scheme than to the dealiased spectral reference.

Appendix figures & tables7 assets

Supplementary material from the paper’s appendix.

Appendix

Panel of Figure 2 — Figure 2: The measurement pipeline. One forward pass per wavenumber returns the empirical symbol $\hat{G}_{\theta}(k)$ of the trained surrogate’s linearized one-step map (Sec. 2.2 ). Multiplying the residual $\hat{G}_{\theta}-G_{\mathrm{exact}}$ by the complex conjugate $\overline{G_{\mathrm{exact}}}$ rotates it into the frame of the exact reference, whose real and imaginary parts are the dissipation and dispersion coordinates (Sec. 2.4 ). Projecting onto the direction $G_{\mathrm{scheme}}-G_{\mathrm{exact}}$ gives the provenance coefficient $\alpha(k)$ (Sec. 2.5 ); differencing the symbols learned from two schemes suppresses the common component of the architecture-intrinsic term $B(k)$ (Sec. 2.6 ).

account	$\alpha_{\mathrm{up}}$	$\alpha_{\mathrm{CN}}$	$\mathrm{Re}\,\Delta_{\mathrm{twin}}$	$\mathrm{Im}\,\Delta_{\mathrm{twin}}$	right angle?
exact-reference-consistent	$0$	$0$	$0$	$0$	no (both at origin)
spectral bias only	$>0$ (spurious)	$\approx 0$	$0$	$0$	no (both radial)
scheme imitation	$1$	$1$	$<0$	$<0$	yes
emulator superiority	$<0$	$<0$	$>0$	$>0$	no (past exact reference)
full-imitation ceiling	$1$	$1$	$-0.241$	$-0.685$	–
observed	$1.00\pm 0.00$	$1.00\pm 0.00$	$-0.241$	$-0.684$	yes

Appendix

Table 4: What each account of surrogate error predicts, and what is observed. Columns are the provenance coefficients of the two surrogates (( 6 )) and the two components of the twin-scheme difference (( 9 )); the last column asks whether the two learned symbols separate along orthogonal directions in the

(\Delta_{\mathrm{diss}},\Delta_{\mathrm{disp}})

plane. The full-imitation ceiling is computed from the analytic symbols on this grid; the observed row is the linear-advection CNN experiment of Sec. 4.1 , mean

\pm

s.d. over 3 seeds. Only scheme imitation predicts a nonzero tangential twin difference.

	radial (dissipation)	% of ceiling	tangential (dispersion)	% of ceiling
CNN (local)	$-0.2406\pm 0.0001$	$100.0\%$	$-0.6841\pm 0.0008$	$99.9\%$
FNO (non-local)	$-0.2402\pm 0.0001$	$99.9\%$	$-0.6837\pm 0.0003$	$99.9\%$
full-imitation ceiling	$-0.2406$	–	$-0.6845$	–

Appendix

Table 5: Twin-scheme difference

\Delta_{\mathrm{twin}}=(\hat{G}_{\theta}^{\mathrm{up}}-\hat{G}_{\theta}^{\mathrm{CN}})\overline{G_{\mathrm{exact}}}

on linear advection, averaged over the upper spectrum, mean

\pm

s.d. over 3 seeds, against the ceiling it would reach under full scheme imitation (

\hat{G}_{\theta}^{(s)}=G_{\mathrm{scheme}}^{(s)}

, at which the learned symbols have no residual

B(k)

). Both architectures saturate both axes; the tangential axis, which a scheme-independent radial

B(k)

cannot produce, is as saturated as the radial one. The two architectures have very different architectural displacements (their exact-trained controls differ by five orders of magnitude in MSE, Table 2 ) yet agree here to

4\times 10^{-4}

, as expected if what survives the difference is the scheme term.

input pool	surrogate	$\alpha$	$D_{\mathrm{scheme}}$	$D_{\mathrm{exact}}$	ratio	artifact
own scheme’s states	upwind	$1.032$ [1.023, 1.040]	0.564	3.02	$5.4\times$	3.48
own scheme’s states	centered	$0.946$ [0.937, 0.954]	0.127	3.40	$26.8\times$	3.46
spectral-reference states	upwind	$2.845$ [2.827, 2.867]	0.754	37.87	$50.2\times$	38.11
spectral-reference states	centered	$2.835$ [2.777, 2.906]	0.720	37.72	$52.4\times$	37.00
mixed (all three)	upwind	$2.836$ [2.782, 2.896]	0.987	3.50	$3.5\times$	3.87
mixed (all three)	centered	$1.968$ [1.869, 2.070]	0.142	3.46	$24.4\times$	3.43

Appendix

Table 6: Forced Burgers, one-step operator with the input distribution held fixed and only the operator varying.

\alpha

is the field-space projection ( 10 ), with

1

meaning the learned operator reproduces its training scheme and

0

meaning it reproduces the dealiased spectral reference; brackets are

95\%

bootstrap confidence intervals over trajectories.

D_{\mathrm{scheme}}

D_{\mathrm{exact}}

, ratio and artifact are as in Table 3 , now computed on the common pool. The first block recomputes the protocol of Table 3 ; the lower blocks evaluate every operator on spectral-reference states and on a mixture of all three solvers’ states. Fixing the input distribution does not remove the preference for the training scheme.

Figure 3: Forced Burgers: operator-level provenance is not visible at the attractor. (a) stationary energy spectra $E(m)$ against mode $m$ (log–log). Solid lines are the three solvers, which separate by scheme as modified-equation theory predicts (centered above the spectral reference above upwind at high $m$ ); dots are the corresponding surrogates run autoregressively to their own steady state, which collapse onto one another and fall below every solver. The shaded band is the fitting range for panel (b). (b) high- $k$ spectral slope fitted over $6\leq m\leq 40$ : the solvers straddle the spectral-reference slope (dotted line), while every surrogate rollout is steeper than every solver, i.e. autoregressive over-smoothing dominates the scheme differences. The scheme signature these same models carry in their one-step operators (Tables 3 and 6 ) is therefore not recoverable from rollout statistics.

	upwind-trained			CN-trained
budget	$\alpha$	MSE exact	MSE scheme	$\alpha$	MSE exact	MSE scheme
0	$1.000\pm 0.000$	$1.1\times 10^{-1}$	$1.4\times 10^{-7}$	$1.000\pm 0.000$	$1.8\times 10^{-1}$	$4.1\times 10^{-8}$
25	$0.008\pm 0.006$	$5.7\times 10^{-3}$	$1.0\times 10^{-1}$	$-0.021\pm 0.010$	$6.6\times 10^{-3}$	$1.8\times 10^{-1}$
50	$0.006\pm 0.001$	$5.6\times 10^{-3}$	$1.0\times 10^{-1}$	$-0.026\pm 0.005$	$6.5\times 10^{-3}$	$1.8\times 10^{-1}$
100	$0.007\pm 0.005$	$5.6\times 10^{-3}$	$1.0\times 10^{-1}$	$-0.027\pm 0.006$	$6.5\times 10^{-3}$	$1.8\times 10^{-1}$
250	$0.013\pm 0.000$	$4.4\times 10^{-3}$	$1.0\times 10^{-1}$	$0.003\pm 0.005$	$4.9\times 10^{-3}$	$1.8\times 10^{-1}$

Appendix

Table 7: Provenance under high-accuracy fine-tuning (mean

\pm

s.d., 3 seeds). Budget is the number of exact-evolution samples used for fine-tuning, against an original scheme corpus of

4000

samples.

\alpha

is measured against the original training scheme:

1

means still scheme-consistent,

0

means the scheme’s error direction is no longer present. MSE exact and MSE scheme are one-step errors against the exact evolution and against the original training scheme; together they show that the model does not interpolate between the two operators but changes target, since error against the exact evolution falls while error against the original scheme rises by six orders of magnitude. Accuracy against the exact evolution saturates at

3.5\times 10^{-3}

, the same floor the exact-trained CNN reaches in Table 2 .

model	$10^{-4}$	$10^{-3}$	$10^{-2}$	$10^{-1}$
cnn/cn	$5.7\times 10^{-5}$	$7.3\times 10^{-6}$	–	$7.1\times 10^{-6}$
cnn/exact	$1.8\times 10^{-4}$	$10.0\times 10^{-5}$	–	$1.0\times 10^{-2}$
cnn/upwind	$3.4\times 10^{-4}$	$5.0\times 10^{-5}$	–	$2.4\times 10^{-3}$
fno/cn	$6.2\times 10^{-5}$	$4.8\times 10^{-6}$	–	$2.5\times 10^{-5}$
fno/exact	$9.1\times 10^{-5}$	$1.2\times 10^{-5}$	–	$7.7\times 10^{-5}$
fno/upwind	$9.7\times 10^{-4}$	$1.2\times 10^{-4}$	–	$7.5\times 10^{-4}$

Appendix

Table 8: Linearity of the symbol probe:

\max_{k}|\hat{G}_{\theta}(\varepsilon)-\hat{G}_{\theta}(10^{-2})|/|\hat{G}_{\theta}(10^{-2})|

over the resolved spectrum, for each architecture/corpus combination. The column at

10^{-2}

is the reference point. The recovered symbol is stable over the tested

10^{-4}

–

10^{-1}

range, with worst-case relative change

1\%

Explore similar work

Sep 21, 2026math.NA

Cost-Accuracy Trade-offs: Neural Operator vs Classical Numerical Solver

Neural operators are data-driven models that learn mappings from inputs that parameterize partial differential equations, such as spatially varying coefficients, initial conditions, forcing terms, boundary conditions, or geometries, to solution fields or quantities of interest. Once trained, they can serve as surrogates for classical numerical solvers in many-query settings that require repeated evaluations for varying inputs. We address the question of when, and then why, neural operator surrogates outperform classical numerical solvers, in terms of cost for a given accuracy. We focus on the post-training, many-query limit, in which data-acquisition and training costs are treated as fixed and fully amortized. Even in this deliberately favorable regime for neural operators, there are regimes in which classical solvers outperform the surrogate models. We compare the cost-accuracy performance of neural operator surrogates and classical numerical solvers through a reproducible benchmark study comparing neural operators with problem-matched classical solvers on representative problems in computational science and engineering, focusing on prediction error, per-query floating-point cost, and wall-clock runtime. Neural operators are most competitive at low-to-moderate accuracy requirements. Their floating-point cost advantage depends strongly on the problem structure, arising when they avoid temporal or nonlinear iterations or predict a reduced quantity of interest rather than a full solution field. Additional wall-clock speedups result from dense tensor operations that are well suited to modern hardware. As the target accuracy is tightened, achieving the required accuracy with neural operators becomes increasingly challenging, and classical solvers outperform surrogates in this regime; thus classical solvers will remain important for verification and high-accuracy computation.

Daniel Zhengyu Huang, Andrew M. Stuart

Beijing International Center for Mathematical Research, Center for Machine Learning Research, Peking University, Beijing, China · California Institute of Technology, Pasadena, CA

Neural Operators Classical Numerical Solver Surrogate Models

Ask PDF

Jun 16, 2026cs.LG

Operator Boosting Produces Pareto-Efficient PDE Surrogates

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 Neural Partial Differential Equation Solvers Surrogate Models Boosting

Ask PDF

May 14, 2026cs.LG

Breakeven complexity: A new perspective on neural partial differential equation solvers

Neural surrogate solvers of partial differential equations (PDEs) promise dramatic speedups over numerical methods, especially in scenarios requiring many solves. However, current accuracy-based evaluations do not fully consider two central issues: (1) neural solvers incur substantial up-front costs for data generation, training, and tuning; and (2) classical solvers can also generate low-fidelity solutions at a sufficiently low simulation cost. To explicitly account for these realities and fully incorporate end-to-end costs, we propose an evaluation framework centered on breakeven complexity, a metric that counts the forward solves before a learned solver is cost-effective relative to an error-equivalent traditional solver. To evaluate this measure, we apply scaling laws to determine how much training budget to allocate to data generation and discuss how to achieve smooth error-matching in diverse settings. We evaluate the breakeven complexity of multiple neural PDE solvers on three PDEs on 2D periodic domains from APEBench and a novel benchmark of flows past multiple obstacles generated by the GPU-native PyFR code. Among other findings, our results suggest that neural PDE solvers become more effective as problems get harder in terms of cost, dimension, rollout, physics regime (e.g. higher Reynolds number), etc.

Yijing Zhang, Nicholas Roberts, Tanya Marwah +1

University of Wisconsin–Madison · Google DeepMind

Neural Partial Differential Equation Solvers Reynolds-Averaged Navier-Stokes Perspective

Ask PDF