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

AI summary

Overview

Research area: Scientific machine learning — specifically neural surrogate models ("neural operators" / neural PDE solvers) that learn to approximate the solutions of partial differential equations from numerical solver data.

Technical level: Advanced. The paper assumes familiarity with discrete Fourier analysis, von Neumann stability/amplification factors, finite-difference discretization error, and neural operator architectures such as Fourier Neural Operators.

One-sentence scope: The paper introduces an empirical "Fourier symbol" diagnostic that measures whether a trained neural PDE surrogate has learned the underlying physics or merely imitated the discretization error of the numerical scheme that produced its training data.

What This Paper Is About

Neural PDE surrogates are usually trained on trajectories generated by a numerical solver, and then evaluated against held-out trajectories from that same solver. The paper argues that this standard protocol is ambiguous: because the test targets are generated by the same scheme as the training targets, a network that has faithfully copied the solver's discretization error (its artificial damping, its phase lag) can score just as well as — or better than — a network that has learned the true physical evolution. The goal is to build a measurement instrument that can tell these two cases apart.

Key Contributions

1. A provenance test built from the empirical Fourier symbol. The authors turn the single-mode probing primitive (described as introduced for operator networks by Gao et al., 2026) into an attribution tool by resolving the learned symbol into dissipation and dispersion coordinates and projecting it onto the axis joining the exact-evolution reference to the training-scheme reference.

2. Twin-scheme differencing. Two identical networks are trained on schemes with orthogonal error signatures — first-order upwind (error almost purely dissipative) and Crank–Nicolson (error purely dispersive) at Courant number σ = 1/2 — with all other experimental factors held fixed. The scheme-independent architectural spectral bias is suppressed in the difference of the two learned symbols, yielding a falsifiable prediction (a specific saturation value) that spectral bias alone cannot reproduce.

3. Application to standard surrogates. The diagnostic is applied to a local periodic five-layer CNN and a non-local FNO on linear advection, at the operator level on forced viscous Burgers with the input distribution held fixed, and the paper shows that provenance and accuracy respond on different scales under high-accuracy fine-tuning.

4. A provenance coefficient α(k). Definition of a scalar projection, α = 1 when the learned symbol coincides with the scheme, α = 0 when it is consistent with the exact reference in the scheme-error direction, and α < 0 for the "emulator superiority" regime in which a surrogate is more accurate than the solver that trained it.

Main Findings

• Linear advection: the learned symbol matches the training scheme, not the exact reference. The upwind-trained CNN reproduces the upwind scheme's amplitude deficit to three decimals (learned |Ĝθ| = 0.395 ± 0.000 against a scheme value of 0.395) while leaving the phase almost untouched (learned lag 0.003 ± 0.003 against a scheme lag of exactly 0). The Crank–Nicolson-trained CNN remains unitary to three decimals and reproduces the scheme's phase error (0.824 ± 0.000 against 0.824). For both, α = 1.000 ± 0.000.

• The exact-trained control reveals architectural displacement. Trained on exact data, the CNN sits near but not at the exact reference (|Ĝθ| = 0.987 ± 0.006, lag 0.025 ± 0.003), which the paper reads as the architecture-intrinsic spectral bias term B(k) made visible.

• The twin-scheme difference saturates its ceiling. Measured over the upper spectrum, Re Δ_twin = −0.2406 ± 0.0001 and Im Δ_twin = −0.6841 ± 0.0008, against a full-imitation ceiling of −0.2406 radial and −0.6845 tangential computed from the analytic symbols. This recovers 100.0% of the available dissipation signal and 99.9% of the available dispersion signal. The abstract and introduction summarize this as the twin-scheme difference reaching more than 99.8% (over 99.8%) of the analytically predicted ceiling.

• The orthogonal geometry is the key signal. In the (Δ_diss, Δ_disp) plane, the upwind surrogate moves along the dissipation axis and the Crank–Nicolson surrogate along the dispersion axis, forming a right angle. The paper states that spectral bias alone predicts both symbols displaced radially, no tangential separation, and Δ_twin = 0 in both components.

• A non-local architecture shows the same behavior. The FNO (retaining all N/2 resolved modes) reproduces upwind's amplitude deficit (0.395 ± 0.000, α = 1.000 ± 0.000) and Crank–Nicolson's phase error (0.825 ± 0.000, α = 1.001 ± 0.000). Its twin difference reaches 99.9% radial and 99.9% tangential of the ceiling, against the CNN's 100.0% and 99.9%.

• The capacity/locality explanation is ruled out in this experiment. Trained on exact data, the CNN cannot fit the non-local exact evolution (MSE 3.3 × 10⁻³) while the FNO fits it to 3.6 × 10⁻⁸ and lands at the exact reference (|Ĝθ| = 0.999 ± 0.000, mean lag 0.001 ± 0.000). An architecture that represents the exact evolution when shown exact data reproduces the scheme when shown solver data.

• Architectural displacements differ wildly but twin differences agree. The CNN and FNO exact-trained controls differ by five orders of magnitude in MSE, yet their twin differences agree to 4 × 10⁻⁴ on both axes (Table 5), which the paper says is what one expects if the quantity surviving the difference is the scheme term.

• Nonlinear Burgers: provenance survives beyond linear dynamics. On forced viscous Burgers, each surrogate's one-step operator action on its own solver's forced steady states is 5–27× closer to its training scheme than to the dealiased spectral reference: upwind has D_scheme = 0.564 vs. D_exact = 3.02 (ratio 5.4×) with a scheme artifact of 3.48; centered has D_scheme = 0.127 vs. D_exact = 3.40 (ratio 26.8×) with a scheme artifact of 3.46.

• Common-input evaluation preserves the effect. With the input distribution held fixed so only the operator varies, the paper reports the learned operator remains substantially closer to its training scheme's action than to the dealiased spectral reference action — 50.2× for upwind (the remaining figures for the other surrogate are cut off in the provided text).

• Autoregressive rollout obscures provenance. Autoregressive Burgers rollout adds a scheme-independent over-smoothing that dominates the attractor-level spectrum, so provenance must be measured at the operator level.

• Provenance and accuracy respond on different scales. Fine-tuning a scheme-trained CNN on exact-evolution samples moves α from 1 to approximately 0 with as few as 25 samples, while exact-reference MSE remains at the architecture-dependent floor.

• The diagnostic is certified as linear. Sweeping ε ∈ [10⁻⁴, 1], the recovered symbol is ε-independent to 1.0 × 10⁻² in relative terms across ε ∈ [10⁻⁴, 10⁻¹] and departs only at ε ~ 1, where the GELU nonlinearity is no longer probed locally. Applied to the numerical solvers themselves, the same procedure returns their analytic symbols to 2 × 10⁻¹⁶.

Methodology in Plain English

The authors start from a simple idea about linear, translation-invariant operations on a periodic grid: such an operation acts on each Fourier mode by multiplying it by a single complex number, called its symbol. This is exactly the same object as the von Neumann amplification factor used in classical numerical analysis to describe how a finite-difference scheme damps or lags each wavenumber.

To measure the symbol a trained network has actually learned, they do the following. For each resolved wavenumber k, they construct a pure cosine wave of small amplitude (ε = 10⁻²) as an input, run it through the network's one-step map once, take a discrete Fourier transform of the output, and divide the k-th output coefficient by the k-th input coefficient. This costs one forward pass per wavenumber, needs no gradients, and requires no access to model internals. For nonlinear problems, the same quantity is read as the same-mode diagonal Fourier response of the network's Jacobian about a base state drawn from the solver's own forced steady states.

Two reference symbols are known in closed form. The exact evolution of linear advection translates each mode without changing its amplitude: G_exact(k) = e^(−iσφ_k). The training scheme applies a different multiplier with a known closed form. The authors use explicit first-order upwind and Crank–Nicolson with centered differences. At the Courant number used (σ = 1/2), upwind's error is entirely a modulus deficit (its phase agrees exactly with the exact evolution) and Crank–Nicolson is unitary with an entirely phase-based error — the two scheme errors separate exactly into a modulus and a phase component.

Because the exact symbol has unit modulus, multiplying the error by its complex conjugate rotates the complex plane so the exact reference sits at 1, and the real and imaginary parts of the rotated residual become a dissipation coordinate (amplitude error) and a dispersion coordinate (phase error). Projecting the learned error onto the direction from exact to scheme gives the provenance coefficient α(k) — a least-squares projection treating ℂ as ℝ².

The confound is that neural networks are architecturally biased toward low frequencies, and that bias damps high wavenumbers, which looks qualitatively like a dissipative scheme's error. A single dissipative-scheme experiment would therefore report α > 0 even for a network that had learned nothing scheme-specific. The fix is twin-scheme differencing: train two identical networks (same architecture, seed sequence, optimizer, schedule, batch size, epoch count, sample count, mesh, time step, and input distribution) on the two schemes whose errors are orthogonal, then subtract their measured symbols. The scheme-independent part of the architectural bias cancels in the difference, leaving a quantity whose full-imitation value can be computed analytically from the two scheme symbols alone. The paper is explicit that target-dependent approximation and optimization effects need not cancel, so this is described as an identification strategy rather than an exact causal decomposition.

The experimental setup is deliberately controlled. For linear advection (u_t + u_x = 0 on [0, 1], N = 64, σ = 0.5), three corpora of 4000 input/target pairs are generated from random band-limited fields with approximately flat spectra: upwind, Crank–Nicolson, and an exact spectral control. Surrogates are a periodic five-layer CNN and an FNO retaining all N/2 resolved modes, trained with Adam and cosine annealing over 400 epochs and three seeds. For forced viscous Burgers (u_t + uu_x = νu_xx + f on [0, 2π], N = 128, ν = 0.04), three method-of-lines RK4 solvers share mesh, time step, viscosity, and forcing and differ only in the convection discretization: dealiased spectral (the high-fidelity reference), first-order upwind, and second-order centered. The surrogates learn the M = 15 substep map.

Why This Matters

Impact on research. The paper makes an argument with consequences for how neural PDE benchmarks are interpreted. Benchmarks such as PDEBench and PDEArena distribute solver-generated trajectories as targets, and standard neural-operator and autoregressive simulator studies evaluate against numerical rollouts. If those rollouts carry the scheme's numerical signature, then low error against them establishes fidelity to the solver, not fidelity to the exact evolution — and fidelity to the exact evolution is only one possible cause of that agreement. The Fourier-symbol diagnostic provides a way to attribute the error rather than merely quantify it, which the paper says the closest prior spectral diagnostic (measuring how far a learned operator sits from the exact operator) does not do: it never compares the exact evolution and the training scheme side by side.

Real-world applications:

• Scientific simulation surrogate development: teams replacing expensive solvers with learned emulators need to know whether the emulator reproduces physics or reproduces the reference solver's artifacts.
• Climate and weather modeling: numerical schemes for advection have well-known damping and phase-lag signatures; a surrogate that inherits them will inherit systematic biases.
• Engineering design and CFD workflows: surrogate models used inside optimization loops can propagate a solver's dispersive phase errors into downstream design decisions.
• Model validation and benchmarking: the diagnostic gives benchmark designers a tool for checking whether their evaluation protocol distinguishes physics fidelity from scheme imitation.

Industry relevance. Any organization substituting a learned surrogate for a classical PDE solver — in aerospace, energy, materials, or geophysics — faces the question this paper formalizes. The measurement requires only forward passes and no model internals, so it can be applied to existing trained models without retraining, and the paper states that provenance and benchmark accuracy are distinct quantities that can respond on different scales. That distinction matters for procurement and validation: a model that scores well on a solver-generated test set may nonetheless be reproducing a scheme's errors rather than the underlying physics.

Future Directions

• Extending beyond the tested scope. The paper states that its evidence covers one-dimensional periodic advection and forced Burgers at modest resolution, two architectures, and three seeds. Whether the same provenance signature appears in higher dimensions, more complex PDEs, or with different solver families is not established here.

• Replacing the numerical "exact" comparator with a true continuum oracle. For Burgers, the paper explicitly notes that the "exact" comparator is a dealiased spectral numerical reference, not a continuum oracle. An open question is how the diagnostic behaves against genuinely exact or analytically known solutions.

• Understanding the incomplete cancellation of B(k). The identification argument requires the architectural spectral-bias term to be scheme-independent, and the paper says it is only approximately so, because the two corpora have different spectral content and therefore different optimization trajectories. A tighter characterization of the residual, target-dependent component remains open.

• Making provenance robust under rollout. The paper finds that autoregressive Burgers rollout adds a scheme-independent over-smoothing that dominates the attractor-level spectrum, so provenance must be measured at the operator level. Whether a rollout-level diagnostic can be constructed is left open.

• Relating fine-tuning sample count to provenance decay. Fine-tuning on 25 exact samples moves α to approximately 0 without removing the architecture-dependent accuracy floor. How α decays with sample count in general, and what governs the accuracy floor, is not resolved in this work.

Target Audience

This paper is most useful to machine learning researchers working on neural operators and physics-informed or data-driven PDE surrogates, and to numerical analysts interested in how classical discretization theory (von Neumann analysis, modified-equation theory) can be turned into a measurement instrument for learned models. Benchmark designers and practitioners who validate surrogate models against solver-generated test data will find the central argument directly applicable to their evaluation protocols. It is less suited to readers without a background in Fourier analysis or numerical discretization, since the diagnostic is defined in those terms throughout.

Authors’ 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.