Neural means and kernel corrections for operator learning

Yitzchak Shmalo

Source-linked AI summary

Neural means and kernel corrections for operator learning

Yitzchak Shmalo

arXiv:2609.00389v1Jul 8, 2026cs.LGmath.PRstat.ML

TL;DR

Operator-learning surrogates must balance neural flexibility against kernel methods’ exact solves and error certificates. This paper combines neural means with exact Matérn corrections of residuals and learned features, then evaluates the pairing on structural mechanics and OCO-2 benchmarks. The hybrid matches the best published structural-mechanics result, improves the OCO-2 emulator on two bands, and identifies when kernel exactness is most useful in feature space.

Problem
The paper examines when neural and kernel operator-learning surrogates should be combined, given their contrasting representations, exactness, and benchmark performance.
Method
The method trains neural network means, stacks multiple architectures, and applies exact Matérn kernel regressions to residuals and learned features.
Results
4.55% structural-mechanics test error matches the best published architecture, while OCO-2 feature-space correction reaches 3.8% versus 40% on the raw state and improves the emulator on two of three bands.
Takeaways & Limitations
Kernel corrections are most valuable either as residual corrections when neural and kernel models tie or inside learned representations when the raw-state kernel trails the network.
Takeaways & Limitations
The distribution-free coverage guarantee is only approximate because calibration uses a split also used for model selection; a disjoint calibration split would make it exact.

Abstract

from arXiv · show

We combine neural network means with exact Matérn kernel regressions of their residuals and of their learned features, and evaluate the pairing on two public emulation problems with published baselines: the structural-mechanics benchmark of de Hoop et al. and the OCO-2 radiative-transfer emulator of Lamminpää et al. On structural mechanics the combination reaches 4.55% test error, matching the best published architecture, and 5.38% against a published 6.49% in the low-data regime. On OCO-2 it improves on the published Gaussian-process emulator on that problem's own test points, outright on two of the three spectral bands; the same kernel that trails the network tenfold on the raw state overtakes it on the network's features, and we measure why (the target's squared native-space norm drops about fortyfold at fixed effective dimension) and prove the mechanism. Where the two families tie instead, the residuals of every architecture we train correlate above 0.86 and their shared component is flat in diversity and sample size, which reads the published plateau as a property of the data. Supporting results include a second-moment identity that predicts stacking outcomes from measured correlations, an optimal-recovery certificate, and a distribution-free coverage band, the only uncertainty signal that survives our tests.

1 Introduction

The paper pairs neural network means with exact Matérn kernel corrections and evaluates the combination across public operator-learning benchmarks. Results show complementary regimes: parity on structural mechanics, representation-level gains on OCO-2, and uncertainty benefits from the supporting theory.

Contribution: The pipeline combines a neural network mean with exact Matérn kernel regressions of residuals and learned features.It is evaluated against published baselines on public emulation problems.
Structural mechanics: 4.55% test error under the 20000-sample structural-mechanics protocol matches PARA-Net, while 5.38% under the 1250-sample protocol compares with a published 6.49%.The high-data result matches the best published architecture within run-to-run noise.
OCO-2: On OCO-2, a feature-space Matérn head reaches 3.8% versus 40% for the same kernel on the raw state.The target’s squared native-space norm falls by about fortyfold at unchanged effective dimension.
OCO-2: The OCO-2 combination improves on the published emulator on both error metrics for two of three spectral bands.On the third band, separate heads win the two metrics.
Theory and uncertainty: Supporting results include a second-moment stacking identity, an optimal-recovery certificate, and distribution-free coverage for the reported uncertainty.The design factor is a poor pointwise error ranking, but split-conformal rescaling is the only uncertainty signal surviving the tests.

2 Benchmark, protocol, and related work

The study uses structural mechanics and OCO-2 as deep benchmarks, supplemented by a broader operator-learning suite and published comparisons. Its protocols preserve task-specific metrics and distinguish published baselines from a single unchanged recipe applied outside the two main problems.

Benchmarks: The structural-mechanics task maps 41-point loads to von Mises stress fields on a regular 41 × 41 grid.The dataset contains 40000 input/output pairs, with the load represented as a vector in R41.
Protocol: Errors are mean relative L2 errors over the test set, computed on grid values in double precision.Trapezoidal quadrature changes reported values by less than 0.02 percentage points.
Protocol: The high-data protocol trains on 20000 samples and tests on the final 20000, while the low-data protocol provides 1250 labels with 1000 used for training.Validation uses 1000 held-out samples in the high-data setting and 250 in the low-data setting.
Published baselines: Published structural-mechanics baselines include PARA-Net at 4.55% in the high-data setting and FNO at 6.49% in the low-data setting.Other comparisons include PCA-Net, DeepONet, and Batlle et al.’s kernel and GP methods.
Related work: Neural operators provide fast, mesh-flexible surrogates, whereas kernel and Gaussian-process methods offer exact solves, convergence theory, and error certificates.The paper positions its hybrid pipeline among neural operators, optimal-recovery kernels, kernel flows, and neural-kernel hybrids.

3 Method

The method builds a neural ensemble on standardized load vectors, stacks its predictions, and applies exact kernel corrections first in input space and then in learned feature space. Validation controls the ensemble and kernel choices, while symmetry augmentation and optional kernel-flow regularization provide additional design components.

Pipeline: The pipeline has three stages: neural means, stacking, and kernel corrections.Networks use the reported per-sample relative L2 error as their training loss.
Neural means: The primary mean is a residual MLP mapping R41 directly to R1681, with penultimate features reused by the feature-space kernel stage.The implementation uses residual SiLU blocks with widths of 1024–1536.
Neural means: A kernel-conditioned refiner receives the kernel prediction alongside the load and learns to correct that prediction rather than reproduce the map from scratch.Out-of-fold kernel predictions keep the training channel independent of each target sample.
Neural means: Additional mean architectures include Fourier neural operators, UNets, and an MSE-trained MLP, while a transformer is excluded after failing to reach the others’ level.The included architectures enter the ensemble and their pairwise residual correlations are measured.
Symmetry: Reflection augmentation and test-time reflection averaging exploit the verified load-to-stress symmetry.The averaging step is guaranteed not to hurt on average and yields small consistent ablation gains.
Stacking: Stacking uses validation-fitted convex weights, with an accepted per-pixel variant that allows grid-dependent weights and an intercept.Weights are initialized from the measured second-moment matrix and polished on held-out data.
Kernel corrections: Kernel ridge regression first corrects stacked residuals with a Matérn-5/2 kernel and then repeats the construction on concatenated penultimate features.Kernel scale and regularization are selected on validation data before refitting on the full correction set.
Scope: The construction is intended to transfer beyond this PDE when low-dimensional input parametrizations and verified symmetries are available.Removing feature correction, stacking, or kernel corrections recovers progressively simpler pipeline variants.

4 Results on the structural-mechanics benchmark

The full pipeline reaches the structural-mechanics benchmark’s published performance ceiling, while residual correlations and scaling tests indicate that the remaining error is largely shared and data-limited. Kernel corrections improve accurate neural means by targeting the load-visible remainder, but input-space uncertainty rankings fail where shared errors dominate.

Main comparison: 4.55% test error under the 20000-sample protocol matches PARA-Net, while the full pipeline’s advantage over other methods exceeds run-to-run variation.The kernel ridge baseline reproduces the published result, supporting protocol alignment; the full pipeline reaches 4.55%, below FNO, PCA-Net, DeepONet, and optimal-recovery kernel results.
Pipeline ablation: 4.55% arises when Matérn regression corrects the neural mean’s residual, whereas stacking alone leaves the single high-data mean unchanged.The neural mean starts at 4.86% and the kernel ridge at 5.19%; residual correction lowers the error to 4.55%.
Architectural diversity and the shared residual: 0.86–0.96 residual correlations across trained networks, and 0.79–0.89 correlations with the kernel predictor, show that architectural diversity does not produce substantially different mistakes.The six-member ensemble’s predicted best-mixture error is 4.97%, while the fitted stack reaches 4.61%, consistent with limited decorrelation.
Architectural diversity and the shared residual: 90% of average residual energy lies in the shared component, which concentrates near loaded-edge corners where the finite-element solution is least reliable.The shared component is larger relative to the local field scale at those corners, linking the residual anatomy to the benchmark’s numerical structure.
Architectural diversity and the shared residual: 4.95%, 4.86%, and 4.86% MLP errors at 3500, 8500, and 19000 samples show a plateau, while kernel correction improves the stack by only a few hundredths at n = 19000.These measurements support interpreting the remaining error as data-limited, though the paper identifies this as an inference without regenerated finite-element data.
Uncertainty quantification: 271× lower fitted squared RKHS norm for residuals than raw stress fields shows that the neural mean leaves a smoother remainder for kernel correction.The corresponding certificate tightens because the bound is linear in the norm, but the input-space design factor does not reliably rank pointwise relative error.
Uncertainty quantification: 91.6% empirical test coverage at a nominal 90% level is the uncertainty result that survives, whereas ensemble disagreement correlates only 0.08 with corrected absolute error.Conformal rescaling provides coverage under exchangeability, but the validation split was also used for selection, so the observed guarantee is empirical rather than exact.

5 The OCO-2 radiative-transfer emulator

On OCO-2, exact Matérn correction is most effective after the neural network learns a useful representation, while metric-specific training and data scaling determine which errors remain.

Representation ladder: 3.82% reduced-metric error: the Matérn kernel on trained features slightly surpasses the network’s 3.99%, while the raw-input kernel scores about 40%.The representation ladder records 40% on raw inputs, 30% after sensitivity rescaling, and 3.82% on learned features.
Metric-specific training: 0.0012 versus 0.0054 radiance RMSE on the leading coefficient: the kernel-flow emulator favors heavily weighted coordinates, while flat training favors trailing coordinates.Training in the radiance metric reallocates accuracy toward the coordinates that dominate the reported radiance error.
Architectural diversity: 3.9%, 16.2%, and 8.0%: independently seeded flat-network ensembles reach same-class floors on O2, WCO2, and SCO2, respectively.The feature-kernel head and final combination move below these floors, indicating that changing the representation matters more than adding same-class seeds.
Band-wise comparison: 3.83% on O2 and 7.96% on SCO2: the final combination beats the kernel-flow emulator on both reduced and radiance metrics.WCO2 is the exception: different heads win its reduced and radiance metrics, so no single combined surrogate is reported.
Data scaling: −0.68 fitted slope: corrected-surrogate error falls from 74% at a few hundred pairs to 4.3% at 17,700 pairs.The raw-state kernel improves much more slowly, making additional training pairs the largest practical lever for this task.
Figure 7: 40% to 3.8%: Figure 7 compares the same Matérn machinery across raw inputs, learned features, and intermediate emulators, with per-coordinate RMSE shown by radiance weight.The central comparison is whether the kernel operates on the original state or on the network’s learned features.

6 Theory

The theory supplies guarantees and diagnostics for each pipeline stage, explaining when stacking helps, how feature representations close kernel gaps, and how residual corrections can be certified. It also identifies benchmark-dependent floors and separates metric from representational advantages.

Symmetry averaging: Proposition 6.1 justifies reflection averaging at test time and reflection-augmented training, with small but strictly positive measured gains for every trained model.The result licenses both symmetry-handling steps used in the pipeline.
Kernel-flow regularization: A kernel-flow regularizer combines empirical fit with a feature-generalization term, training the neural mean to produce features suited to the subsequent kernel correction.The empirical-risk term measures fit, while β e2 measures generalization of a kernel method on learned features.
Stacking and ensemble floors: No number of additional members with the same error level and mutual correlations moves ensemble root-mean-square relative error below ē√ϱ.The second-moment decomposition predicts this floor from member errors and pairwise residual correlations.
Stacking and ensemble floors: On structural mechanics, residual correlations of 0.86–0.96 and flat shared error explain why the measured stack reaches a benchmark floor near the best single member.The evidence supports reading the plateau as a property of the finite-element data rather than a surrogate class.
Metric and representation: The measured factor of 42 in native-space norm, at nearly unchanged effective dimension, quantifies the representational advantage of learned features.The pullback interpretation says the raw kernel carries the warped map, whereas the feature kernel carries a simpler factorized target.
Residual correction and uncertainty: The optimal-recovery bound factors corrected-surrogate error into the residual native-space norm and a kernel-design factor Pλ, interpreted as posterior standard deviation.The bound is pathwise and valid for every fixed mean, dataset, and input, but Pλ ranks actual error poorly on the benchmark.

7 Discussion

The discussion frames the method as adaptive across two regimes: coupling helps when neural and kernel models tie, while feature-space correction helps when the kernel trails on raw inputs. Residual correlations, error ratios, and benchmark-specific limitations determine how far performance can move.

Relation to neural-mean Gaussian processes: The correction stage differs from jointly trained neural-mean Gaussian processes by fitting the neural mean first, validating the kernel separately, using an exact solve, and adding feature-space correction.The measured drop in fitted RKHS norm is presented as the quantitative rationale for correcting accurate means.
Limitations: The exact kernel solve relies on one-dimensional input functions and fixed geometry, while the certified bound depends on an empirically lower-bounded residual RKHS norm.Split-conformal calibration is described as the assumption-light complement, and reported error levels remain specific to the dataset and protocol.
Decision diagnostics: Pairwise residual correlations set the ensembling floor, while the member error ratio supplies a keep-or-drop decision for candidate models.These diagnostics are presented as applicable beyond the two evaluated datasets.
Outlook: Refined meshes, higher-order elements, or output grids resolving corner singularities could move the measured structural-mechanics error floor.The discussion identifies dataset regeneration as a route to testing surrogate quality rather than the benchmark’s own discretization limits.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process

The author reports substantial use of generative AI tools during manuscript preparation and retains responsibility for the manuscript’s content.

Declaration: Generative AI assisted software development, dataset preparation, emulator-code adaptation, and experiment implementation.The author states that the resulting work was reviewed and that responsibility for the manuscript remains with the author.

Data and code availability

The paper provides code, experiment summaries, configurations, and data-access information for auditing the reported results.

Availability: The repository contains code, per-run summaries, and exact experiment configurations, while benchmark and OCO-2 resources are identified through their accompanying data records and OSF project.The listed resources cover structural-mechanics, Helmholtz, Navier–Stokes, advection, OCO-2, and ClimSim materials.

A Proofs

The proofs establish symmetry-based averaging, interpolation and expectation identities, ensemble-risk bounds, simplex discretization guarantees, and optimal weighting results. Together, they formalize the pipeline’s averaging, correction, and stacking components.

Symmetry: Equivariance under an involution and invariant loss imply the reflected prediction relation used by reflection averaging.The proof applies equivariance twice and uses invariance of the input distribution for expectation-based conclusions.
Kernel interpolation: Strict positive definiteness makes the restricted-kernel interpolant unique and exact on the selected index set.The interpolation terms indexed by the selected set vanish, establishing the first claim used in the proof.
Expected error: The proof exchanges differentiation and expectation under local integrable Lipschitz conditions while analyzing the held-out kernel error.Conditional counting identifies the number of terms in the relevant sum.
Ensemble bounds: Weighted ensemble error is bounded by the weighted member errors, with the risk inequality following after normalization and expectation.A uniform per-sample loss bound also yields a bounded loss for the combined predictor.
Simplex generalization: Uniform risk control over the simplex is obtained by Lipschitzness, grid discretization, Hoeffding concentration, and a union bound.The grid cardinality determines the finite-complexity term in the resulting guarantee.
Simplex generalization: An empirical simplex minimizer is compared with the population minimizer through twice the uniform empirical-to-population deviation.The proof chooses k = m(M −1) to control the discretization complexity.
Optimal stacking: For two members, the optimal mixture is interior precisely when residual correlation is below the member-error ratio; otherwise the second member is dropped.For exchangeable members, the covariance structure yields the corresponding multi-member result.
Coordinatewise weighting: With diagonal metrics, coordinatewise minimization gives independent optimal weights, and shared weights are dominated unless they already match those coordinatewise optima.The shared-weight case is a restricted subset of the feasible coordinatewise problems.

A.7 Proof of Proposition 6.8

The proof establishes pullback-kernel structure and derives approximation-rate comparisons for isotropic versus adapted Matérn kernels, alongside a residual-based uncertainty bound valid beyond interpolation.

Pullback RKHS: The composition operator Tf=f∘φ has bounded point evaluations because (Tf)(x)=⟨f,k(·,φ(x))⟩Hk.This identifies evaluation on the composed space with the original RKHS representer.
Pullback RKHS: The image space Hkφ is an RKHS with kernel kφ(x,x′)=k(φ(x),φ(x′)) and quotient norm bounded by the original Hk norm.For h∈Hk, the minimum-norm characterization gives ∥h∘φ∥Hkφ≤∥h∥Hk.
Approximation rates: Matérn kernel ridge regression obeys a native-space error estimate controlled by fill distance, with effective dimension determining the sampling rate up to logarithmic factors.For n samples in effective dimension d′, the fill distance scales as n^-1/d′ up to the absorbed logarithmic factor.
Approximation rates: For F=g(A·), isotropic sampling uses n^-1/d, whereas an adapted kernel can achieve n^-1/r when trailing singular directions fall below the resolution threshold.If A is close to full rank, both adapted and isotropic fill distances scale as n^-1/d and the rates coincide.
Uncertainty bound: The residual componentwise bound holds for every λ≥0, and its factor ePλ is no larger than the reported posterior standard deviation Pλ.Thus smoothing can tighten this particular bound relative to the posterior standard deviation.

A.10 Proof of Proposition 6.10

The proof uses exchangeability of normalized validation and test scores to derive finite-sample coverage bounds for the conformal prediction band.

Exchangeability: The normalized scores si=∥e(ũi)∥2/Pλ(ũi) and s⋆=∥e(u⋆)∥2/Pλ(u⋆) are exchangeable when inputs are exchangeable and independent of the fitted quantities.The prediction event is equivalent to s⋆≤q, where q is the selected validation order statistic.
Coverage guarantee: Uniform rank of the test score among m+1 augmented scores yields the lower coverage bound k/(m+1), with k=⌈(1−α)(m+1)⌉.The event s⋆≤q contains the event that the test score has one of the k smallest ranks.
Coverage guarantee: With almost surely distinct scores, coverage is also bounded above by 1−α+1/(m+1).Without ties, the threshold event and the rank event coincide exactly.

A.11 Proof of Corollary 6.12

The corollary derives the effective-dimension limit for kernel degrees of freedom by transferring the spectrum of the normalized Gram matrix to the sample covariance matrix and applying Marchenko–Pastur convergence.

Spectral reduction: The nonzero eigenvalues of K/n=XX⊤/n equal those of S=X⊤X/n, so zero eigenvalues do not affect deff(λ).This reduces the degrees-of-freedom calculation to the p-dimensional sample covariance spectrum.
Asymptotic limit: The empirical spectral distribution of S converges almost surely to the Marchenko–Pastur law with ratio γ and scale σ².Because (x+λ)^-1 is bounded and continuous on [0,∞), the corresponding Stieltjes transforms converge.
Asymptotic limit: Evaluating the Marchenko–Pastur self-consistent equation requires selecting the root whose Stieltjes transform is positive.The sign choice is fixed by positivity of m(−λ).

B Implementation details

The implementation combines neural architectures, kernel residual and feature stages, stacking procedures, and conformal uncertainty calibration under a fixed single-laptop experimental setup.

Experimental setup: All experiments use one laptop with an NVIDIA RTX PRO 2000, 16 CPU cores, 64 GB RAM, and fixed structural-mechanics train/test blocks of 20,000 samples each.Network training uses single precision, while kernel solves and error computations use double precision.
Neural architectures: The reported neural means include FNO, UNet, MLP, and MSE-variant architectures, while the Transformer is excluded because available hardware prevented comparable training.The Transformer is described as a follow-up architecture for stable hardware rather than a reported ensemble member.
Neural architectures: The Refiner concatenates the load with the kernel method’s stress prediction and uses four-fold out-of-fold kernel predictions during training.Reflection augmentation consistently permutes the kernel channel and target by the row-reversal index.
Stacking: Global stacking uses convex weights initialized from the measured second-moment matrix, whereas per-pixel affine weights use ridge regularization and held-out validation acceptance.Per-pixel weights are refitted on the full validation split only after outperforming global weights on the held-out half.
Kernel stages: Matérn-5/2 kernel stages tune scale and nugget grids on validation, refit selected settings on full training data, and apply the feature stage to penultimate activations.The pure-kernel baseline uses a separate scale and nugget grid directly at n=19000.
Uncertainty: Uncertainty uses posterior standard deviation from Cholesky solves, then scales it by the validation 0.9-quantile of normalized error norms.Coverage is measured as the fraction of test samples whose absolute error norm lies below the scaled band.