Neural means and kernel corrections for operator learning

# Neural means and kernel corrections for operator learning
Source: [https://arxiv.org/html/2609.00389](https://arxiv.org/html/2609.00389)
Yitzchak ShmaloEinstein Institute of Mathematics, The Hebrew University of Jerusalem, Jerusalem, Israel\.yitzchak\.shmalo@gmail\.com\. Code, data pointers, and the per\-run summaries behind every number reported here are at[https://github\.com/yspennstate/neural\-means\-kernel\-corrections](https://github.com/yspennstate/neural-means-kernel-corrections)\.

\(July 2026\)

###### Abstract

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\.

## 1Introduction

Operator learning constructs fast surrogates for the solution maps of parametric partial differential equations and physical forward models\. Neural operators\[[8](https://arxiv.org/html/2609.00389#bib.bib6),[12](https://arxiv.org/html/2609.00389#bib.bib4),[11](https://arxiv.org/html/2609.00389#bib.bib5)\]learn their own representations; kernel and Gaussian\-process methods\[[2](https://arxiv.org/html/2609.00389#bib.bib1),[15](https://arxiv.org/html/2609.00389#bib.bib11),[16](https://arxiv.org/html/2609.00389#bib.bib9)\]come with exact solves and error certificates and a strong record on the same benchmarks\. This paper combines the two: a neural network mean, trained in the metric the application reports, with an exact Matérn kernel regression of its residual and of its learned features\. We evaluate the combination in depth on two public problems with published baselines — the structural\-mechanics benchmark ofde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]and the OCO\-2 radiative\-transfer emulation problem ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\]— and, with one unchanged recipe, across the seven\-benchmark operator\-learning suite assembled byBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\]from the problems ofde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]andLuet al\.\[[13](https://arxiv.org/html/2609.00389#bib.bib3)\]\. Table[1](https://arxiv.org/html/2609.00389#S1.T1)collects the published numbers and ours\.

Table 1:The two problems studied in depth \(OCO\-2, structural mechanics\) and the rest of the operator\-learning suite, with the best published mean relativeL2L^\{2\}test error and ours\. OCO\-2 rows compare against the Gaussian\-process emulator ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\]on its own test points, in that paper’s reduced\-coefficient metric; Section[5](https://arxiv.org/html/2609.00389#S5)reports both of its metrics, on which our O2 and strong\-CO2surrogates win both while the weak\-CO2metrics are won by two different heads\. The published suite numbers are per\-problem methods tuned by their authors; ours outside the two problems studied in depth come from one recipe run unchanged \(a Matérn kernel at the median length scale, fit capped at 6000 points, a residual MLP mean, outputs reduced by PCA, the stage selected on validation\)\. Burgers and Darcy use slightly smaller training splits than their baselines \(800 and 824 pairs\) and Darcy a different grid, from the copies we could obtain\.The two problems studied in depth sit at opposite ends of the neural\-versus\-kernel spectrum, and the pair is the point\. On structural mechanics the families tie: the pipeline reaches 4\.55% under the 20000\-sample protocol, matching the best published architecture \(PARA\-Net, 4\.55%\) within run\-to\-run noise, and 5\.38% under the 1250\-sample protocol against a published best of 6\.49%\. What the coupling adds there is mostly understanding\. The residuals of every architecture we train correlate between 0\.86 and 0\.96, their shared component concentrates where the finite element data are least reliable and is flat in ensemble diversity and training\-set size, and so the evidence reads the published plateau near 4\.5% as a property of the benchmark’s data rather than of any surrogate class \(Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)\)\. Absent regenerated finite\-element data this remains an inference, but it is the reading every measurement we made supports\. On the rest of the suite the same recipe sorts the problems by the same logic: where the map is smooth and data are plentiful \(Helmholtz, Navier–Stokes, Darcy\) the kernel stage wins the pipeline’s internal selection and lands within a factor of two of the tuned per\-problem kernels ofBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\], and on Burgers the corrected mean sits near the Fourier neural operator\.

On OCO\-2 the same components separate by an order of magnitude, and the pipeline’s job changes: not to couple two comparable models but to move the kernel’s exactness onto the network’s representation\. Per spectral band of the satellite, the task maps a reduced atmospheric state to a reduced radiance spectrum, and the baseline is the kernel\-flow emulator ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\], whose stored predictions we score on the same test points\. An exact Matérn head on the trained network’s features reaches3\.8%3\.8\\%where the same kernel on the raw state reaches40%40\\%; the effective dimension of the two Gram matrices is unchanged, the target’s squared native\-space norm falls by a factor of about forty, and Proposition[6\.8](https://arxiv.org/html/2609.00389#S6.Thmtheorem8)gives the mechanism\. Trained in the reported metric and combined per output coordinate, the result improves on the emulator on both of that problem’s error metrics at once on two of the three bands; on the third each metric is won by a separate head \(Section[5](https://arxiv.org/html/2609.00389#S5)is precise about which\)\.

Each stage of the pipeline carries one supporting result: a symmetrization lemma behind the reflection averaging, a second\-moment identity that predicts stacking outcomes from measured residual correlations, an optimal\-recovery bound‖G−m‖K​Pλ​\(u\)\\\|G\-m\\\|\_\{K\}\\,P\_\{\\lambda\}\(u\)for the corrected surrogate, a pullback identity and a rate for the feature kernel, and a distribution\-free coverage statement for the reported uncertainty\. The last one matters because the design factorPλP\_\{\\lambda\}, though a valid bound, turns out to be a poor pointwise ranking of the actual error, and its split\-conformal rescaling is the only uncertainty signal that survives our tests \(Section[4\.4](https://arxiv.org/html/2609.00389#S4.SS4)\)\.

Section[2](https://arxiv.org/html/2609.00389#S2)describes the problems, protocols, and published results\. Section[3](https://arxiv.org/html/2609.00389#S3)specifies the pipeline\. Sections[4](https://arxiv.org/html/2609.00389#S4)and[5](https://arxiv.org/html/2609.00389#S5)report the two studies, including the data\-scaling experiments\. Section[6](https://arxiv.org/html/2609.00389#S6)states the supporting theory, one result per stage, with proofs in Appendix[A](https://arxiv.org/html/2609.00389#A1), and Section[7](https://arxiv.org/html/2609.00389#S7)discusses limitations\.

## 2Benchmark, protocol, and related work

### 2\.1Problem and data

The dataset originates in the cost\-accuracy study ofde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]and is the one distributed withBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\]andMoraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\]\. An isotropic elastic plate occupiesD=\(0,1\)2D=\(0,1\)^\{2\}; the displacement fieldwwsatisfies∇⋅σ=0\\nabla\\cdot\\sigma=0with a fixed constitutive law, displacement conditions on the bottom and lateral parts of the boundary, and a prescribed normal tractiont¯\\bar\{t\}on the top edgeΓt=\[0,1\]×\{1\}\\Gamma\_\{t\}=\[0,1\]\\times\\\{1\\\}\. The quantity of interest is the von Mises stress fieldσv\\sigma\_\{v\}onDD\. The learning task is the mapG:t¯↦σvG:\\bar\{t\}\\mapsto\\sigma\_\{v\}\. Loads are drawn from the Gaussian field𝒢​𝒫​\(100,4002​\(−Δ\+32​I\)−1\)\\mathcal\{GP\}\\big\(100,\\,400^\{2\}\(\-\\Delta\+3^\{2\}I\)^\{\-1\}\\big\)with homogeneous Neumann boundary conditions for the Laplacian; outputs are finite element solutions interpolated to a regular41×4141\\times 41grid, and the load is sampled at4141points\. The distributed file contains 40000 input/output pairs; the input array stores the load broadcast along the second grid coordinate, which we verified is an exact copy \(maximal deviation0across all 40000 samples\), so all our methods consume the load as a vector inℝ41\\mathbb\{R\}^\{41\}\.

Errors are mean relativeL2L^\{2\}errors over the test set,

1Ntest​∑n‖G^​\(un\)−G​\(un\)‖2‖G​\(un\)‖2,\\frac\{1\}\{N\_\{\\mathrm\{test\}\}\}\\sum\_\{n\}\\frac\{\\\|\\widehat\{G\}\(u\_\{n\}\)\-G\(u\_\{n\}\)\\\|\_\{2\}\}\{\\\|G\(u\_\{n\}\)\\\|\_\{2\}\},computed on grid values in double precision\. Quadrature weighting \(trapezoidal instead of plain vector norms\) changes our reported numbers by less than 0\.02 percentage points, and we report the plain\-norm convention of the prior work\.

### 2\.2Protocols and published results

Two regimes appear in the literature, and we follow both\. In the*high\-data*protocol the first 20000 samples are available for training and the last 20000 form the test set;de Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]report, at 20000 training samples, 4\.55% \(PARA\-Net\), 4\.67% \(PCA\-Net\), 4\.76% \(FNO\) and 5\.20% \(DeepONet\), andBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\]report 5\.18% for their optimal\-recovery kernel method \(Matérn/rational quadratic; 27\.11% for the linear kernel\) on the same task\. Within the 20000\-sample budget we hold out the last 1000 samples \(of a fixed permutation\) for model selection and stacking, and train on the remaining 19000, so no method of ours sees more than the 20000 training samples used by the baselines\. In the*low\-data*protocol ofMoraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\], 1250 samples are available and the same 20000\-sample test set is used; their table gives 8\.70% \(DeepONet\), 6\.62% \(FNO\), 6\.95% \(the kernel method ofBatlleet al\.[2](https://arxiv.org/html/2609.00389#bib.bib1)\), 6\.74% \(their zero\-mean GP\), 7\.12% \(DeepONet\-mean GP\) and 6\.49% \(FNO\-mean GP\)\. In this regime we carve the validation split \(250 samples\) out of the 1250, so training uses 1000 samples and every choice made by the pipeline is informed by the 1250 available labels only\.

### 2\.3A data\-driven check of the mirror symmetry

The continuous problem is invariant underx1↦1−x1x\_\{1\}\\mapsto 1\-x\_\{1\}: the domain and boundary partition are symmetric, and the input law has a symmetric covariance\. On the grid, reflecting the load \(SS\) should reflect the stress field along the first grid coordinate \(TT\), i\.e\.G​\(S​u\)=T​G​\(u\)G\(Su\)=TG\(u\)\. Rather than assume this, we test it on the data\. For each of the first 200 samples we searched all 40000 loads for the nearest neighbor of the reflected loadS​uiSu\_\{i\}and compared the corresponding outputs\. For the five best\-matching pairs, the input mismatch‖S​ui−uj‖/‖uj‖\\\|Su\_\{i\}\-u\_\{j\}\\\|/\\\|u\_\{j\}\\\|ranges over0\.270\.27–0\.320\.32, while the output mismatch‖T​vi−vj‖/‖vj‖\\\|Tv\_\{i\}\-v\_\{j\}\\\|/\\\|v\_\{j\}\\\|ranges over0\.0850\.085–0\.140\.14; reflecting along the second coordinate instead gives output mismatches of order one\. Outputs of near\-mirror inputs are far closer than the inputs themselves, which is what equivariance plus a Lipschitz solution map predicts, and the effect singles out the correct output reflection axis\. A complementary check on the input law: the empirical mean and covariance of the training loads are invariant under reflection to within 1\.1% and 1\.5% respectively, as required for Proposition[6\.1](https://arxiv.org/html/2609.00389#S6.Thmtheorem1)\. Consistently with all of this, reflection averaging at test time improves every model we trained \(Section[4](https://arxiv.org/html/2609.00389#S4)\), as Proposition[6\.1](https://arxiv.org/html/2609.00389#S6.Thmtheorem1)says it must on average when the symmetry holds\.

### 2\.4Related work

The benchmark sits at the meeting point of three lines of work\.*Neural operators*learn maps between function spaces: DeepONet\[[12](https://arxiv.org/html/2609.00389#bib.bib4)\]with its branch/trunk factorization, the Fourier neural operator\[[11](https://arxiv.org/html/2609.00389#bib.bib5)\]and the broader neural\-operator framework\[[8](https://arxiv.org/html/2609.00389#bib.bib6)\], and the reduced\-basis PCA\-Net and PARA\-Net architectures assembled for the cost–accuracy study ofde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]\. These methods are fast and mesh\-flexible but do not by themselves come with error control, and it is their numbers on this problem that define the plateau we start from\.*Kernel and Gaussian process methods*for operator learning approach the same maps through reproducing kernels: the optimal\-recovery framework ofBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\], with its convergence theory and a\-priori bounds, is the most direct comparison, and it is competitive with neural operators on most of the benchmarks it considers\. Data\-adapted kernels learned by cross\-validation\[[4](https://arxiv.org/html/2609.00389#bib.bib34)\]and vector\-valued kernel formulations extend the reach of this family\.*Hybrids*that place a neural network and a kernel in the same estimator are the closest relatives of our pipeline:Moraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\]use a neural operator as the mean of a Gaussian process and fit the two together, and this is the work we most directly build on, differing in that we fit the mean first and the kernel after, tune by cross\-validation rather than marginal likelihood, and solve the correction exactly at nineteen thousand points\.

Two further threads inform the design\. The kernel\-flow method\[[17](https://arxiv.org/html/2609.00389#bib.bib7)\]learns a kernel from data by minimizing a half\-sample cross\-validation loss, and its use as a regularizer for the inner layers of a network\[[23](https://arxiv.org/html/2609.00389#bib.bib8)\]is exactly the term we analyze in Proposition[6\.2](https://arxiv.org/html/2609.00389#S6.Thmtheorem2); kernel flows have since been used at scale as emulators for physical forward models, for instance in atmospheric radiative transfer retrievals\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\]and in the inference of convective\-storm structure from satellite observations\[[18](https://arxiv.org/html/2609.00389#bib.bib23)\], where presenting the input in its physical parametrization and adapting the kernel to data are what make the emulator accurate\. Our pipeline follows the same instinct, with the smoothing of the target performed by a neural ensemble rather than by a hand\-chosen reduction\. Finally, the effective\-dimension reading of the kernel solve \(Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)\) draws on the random\-matrix description of kernel Gram spectra\[[7](https://arxiv.org/html/2609.00389#bib.bib20)\]and on the classical Marčenko–Pastur\[[14](https://arxiv.org/html/2609.00389#bib.bib14)\]and spiked\-covariance\[[1](https://arxiv.org/html/2609.00389#bib.bib15)\]laws; the same effective dimension governs the generalization of kernel ridge regression\[[3](https://arxiv.org/html/2609.00389#bib.bib19)\]\.

## 3Method

The pipeline has three stages: neural means, stacking, and kernel corrections\. All components operate on the load vectoru∈ℝ41u\\in\\mathbb\{R\}^\{41\}\(standardized by training statistics\) and produce stress fields inℝ41×41\\mathbb\{R\}^\{41\\times 41\}; all networks are trained with the reported metric as the loss, i\.e\. the per\-sample relativeL2L^\{2\}error in original units, which we found mildly but consistently better than mean squared error on normalized targets\.

### 3\.1Neural means

*Residual MLP\.*The primary mean is a residual multilayer perceptron: three or five residual SiLU blocks of width 1024–1536 mappingℝ41→ℝ1681\\mathbb\{R\}^\{41\}\\to\\mathbb\{R\}^\{1681\}directly\. This is essentially PARA\-Net\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]with a modern training recipe, and it also supplies the penultimate features used by the feature\-space kernel stage\.

*Kernel\-conditioned refiner\.*The second mean is the same residual network given an extra input channel: the kernel method’s predicted stress field for the same load, concatenated with the load\. It therefore learns to correct the kernel predictor rather than to reproduce the map from scratch, and it is the strongest single model in our study\. During training the refiner reads*out\-of\-fold*kernel predictions \(four\-fold, so the kernel channel is never fit on the target sample\), and at test time the full\-data kernel prediction; this keeps the channel honest\.

*Other instances\.*The mean slot is not tied to a particular architecture\. We also use a Fourier neural operator\[[11](https://arxiv.org/html/2609.00389#bib.bib5)\]that consumes the broadcast load as a field, a UNet on the same field representation, and an MSE\-trained variant of the MLP; all are drop\-in means, all enter the ensemble of Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3), and their pairwise error correlations are one of the measurements of this paper\. We additionally implement an encoder–decoder transformer\[[19](https://arxiv.org/html/2609.00389#bib.bib17),[6](https://arxiv.org/html/2609.00389#bib.bib18)\]that tokenizes the 41 load samples and decodes the field by cross\-attention at the grid points; it is a natural fourth architecture family, but we were unable to train it to the level of the others on the hardware available for this study and report the ensemble without it \(Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)returns to this\)\.

All means are trained with AdamW and a cosine schedule, batches of 128–256, for up to a few hundred epochs \(Appendix[B](https://arxiv.org/html/2609.00389#A2)lists every hyperparameter\)\. During training each batch is reflected with probability1/21/2\(load reversed, target field reflected along the first grid coordinate\); at test time each model is evaluated as12​\(f​\(u\)\+T​f​\(S​u\)\)\\frac\{1\}\{2\}\(f\(u\)\+Tf\(Su\)\)\. Proposition[6\.1](https://arxiv.org/html/2609.00389#S6.Thmtheorem1)guarantees the test\-time step cannot hurt on average, and the ablation confirms small consistent gains from both steps\. Optionally, the training loss carries a kernel\-flow termβ​e2\\beta\\,e\_\{2\}computed from the batch \(Section[6\.2](https://arxiv.org/html/2609.00389#S6.SS2)\) with a Gaussian kernel on pooled penultimate features whose log\-bandwidth is trained jointly; this variant appears in the ablation\.

### 3\.2Stacking

WithMMtrained modelsf1,…,fMf\_\{1\},\\dots,f\_\{M\}we form the convex combinationm​\(u\)=∑mwm​fm​\(u\)m\(u\)=\\sum\_\{m\}w\_\{m\}f\_\{m\}\(u\),w∈ΔM−1w\\in\\Delta^\{M\-1\}, minimizing the validation relative error; the weights are initialized at the simplex minimizer of the measured second\-moment matrix of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)and polished by a short search on the 1000 held\-out samples \(250 in the low\-data protocol\)\. A per\-pixel variant allows the weights \(plus an intercept\) to vary over the grid, ridge\-fitted on half the validation split and accepted only when it beats the global weights on the other half; it is the variant the final pipeline uses\. Stacking on a held\-out split rather than uniform averaging costs nothing, guards against a weak member\[[22](https://arxiv.org/html/2609.00389#bib.bib30)\], and is statistically almost free at this size \(Proposition[6\.3](https://arxiv.org/html/2609.00389#S6.Thmtheorem3)\)\.

### 3\.3Kernel corrections

The stacked mean is then corrected by kernel ridge regression of its residuals\. LetX=\(u1,…,un\)X=\(u\_\{1\},\\dots,u\_\{n\}\)be the training loads \(standardized coordinatewise\),R∈ℝn×1681R\\in\\mathbb\{R\}^\{n\\times 1681\}the matrix of training residualsvi−m​\(ui\)v\_\{i\}\-m\(u\_\{i\}\), andkka Matérn\-5/25/2kernelk​\(u,u′\)=κ​\(‖u−u′‖/s\)k\(u,u^\{\\prime\}\)=\\kappa\\big\(\\\|u\-u^\{\\prime\}\\\|/s\\big\)\. The correction is

r^​\(u\)=k​\(u,X\)​α,α=\(K\+n​λ​I\)−1​R,\\widehat\{r\}\(u\)=k\(u,X\)\\,\\alpha,\\qquad\\alpha=\(K\+n\\lambda I\)^\{\-1\}R,with\(s,λ\)\(s,\\lambda\)chosen on the validation split from a small grid \(ssas a multiple of the median pairwise distance,λ∈\{10−7,10−5,10−3\}\\lambda\\in\\\{10^\{\-7\},10^\{\-5\},10^\{\-3\}\\\}; the grid is searched on an 8000\-sample subsample and the chosen pair is refit on all ofXX\)\. The corrected surrogate ism\+r^m\+\\widehat\{r\}\. A second corrective stage repeats the construction with the kernel evaluated on deep features, the concatenated penultimate activations of the ensemble members \(reflection\-averaged, standardized\), fitted to the residuals ofm\+r^m\+\\widehat\{r\}; it contributes a smaller improvement and is included when it helps on validation\. Atn=19000n=19000the exact solve is a single Cholesky factorization of a19000×1900019000\\times 19000matrix, under half a minute in double precision on a laptop CPU \(measured in Section[4\.6](https://arxiv.org/html/2609.00389#S4.SS6)\), and prediction is two matrix products; no inducing points, preconditioners, or stochastic solvers are involved\. The Gaussian process reading of the same formulas supplies the posterior standard deviationPλ​\(u\)P\_\{\\lambda\}\(u\)of \([2](https://arxiv.org/html/2609.00389#S6.E2)\) at negligible extra cost \(one triangular solve per evaluation batch\), which is the uncertainty estimate studied in Section[4](https://arxiv.org/html/2609.00389#S4)and certified by Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)\.

Two remarks on scope\. Nothing in the construction uses properties of this particular PDE beyond the verified reflection symmetry, so the recipe \(native input parametrization, symmetrized accurate means, validated kernel correction of the residual, posterior standard deviation as certificate\) transfers to other operator learning problems with low\-dimensional input parametrizations\. And the pipeline degrades gracefully: dropping the feature stage, the stacking, or the corrections recovers progressively simpler methods whose individual numbers appear in the ablation table\.

## 4Results on the structural\-mechanics benchmark

![Refer to caption](https://arxiv.org/html/2609.00389v1/x1.png)Figure 1:Example predictions of the full pipeline\. Top: a median\-error test case; bottom: a case at the 98th error percentile\. Columns: input load, true von Mises stress, prediction, and pointwise absolute error \(note the smaller scale\)\. Errors concentrate in the high\-traction corner, as also reported byMoraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\]\.### 4\.1Main comparison

Table[2](https://arxiv.org/html/2609.00389#S4.T2)reports the high\-data protocol\. The published numbers are quoted fromde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]\(as tabulated byBatlleet al\.[2](https://arxiv.org/html/2609.00389#bib.bib1)\) and fromBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\]; our rows are computed under the split of Section[2](https://arxiv.org/html/2609.00389#S2), which gives our methods strictly no more training information than the baselines had\. The kernel ridge row reproduces the published kernel result almost exactly \(5\.19% against 5\.18%\), which we take as evidence that the protocols are aligned; the gap opened by the rest of the pipeline is therefore not an artifact of evaluation choices\. The full pipeline reaches 4\.55%\. This is below every published method other than PARA\-Net \(4\.55%\), which it matches at the reported two\-decimal precision; since our own stacking configurations vary by a hundredth or two among themselves \(the ablation below moves between 4\.65% and 4\.67% under changes that should not matter\), we read the two as tied and do not claim to separate them\. The margin over the FNO \(4\.76%\), PCA\-Net \(4\.67%\), DeepONet \(5\.20%\) and the optimal\-recovery kernel \(5\.18%\) is several times that spread\. Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)argues the tie is better read as a floor the data impose than as a coincidence of tuning\.

Table 2:Structural mechanics, high\-data protocol \(20000 training samples, 20000 test\)\. Mean relativeL2L^\{2\}test error\. Rows above the rule are published results on the same data; rows below are ours \(mean over the fixed split; TTA denotes reflection averaging at test time\)\.MethodError \(%\)ParametersDeepONet\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]5\.20–FNO\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]4\.76–PCA\-Net\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]4\.67–PARA\-Net\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]4\.55–Optimal\-recovery kernel\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\]5\.18–Kernel ridge on loads \(ours\)5\.19–Residual MLP \+ reflection TTA4\.864\.9M\+ residual kernel correction \(full pipeline\)4\.554\.9MTable 3:Low\-data protocol \(1250 training samples, 20000 test\)\. Published rows fromMoraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\]\. Our pipeline sees only the 1250 labels \(250 of them used as the validation split\)\.
### 4\.2Ablation

Table[4](https://arxiv.org/html/2609.00389#S4.T4)builds the pipeline up one component at a time under the high\-data protocol\. The kernel ridge on loads and the reflection\-averaged neural mean start at 5\.19% and 4\.86% respectively\. Stacking is trivial with a single high\-data mean and leaves it unchanged; the value of the combination appears in the correction stage, where regressing the neural mean’s residual with the validated Matérn kernel lowers the error to 4\.55%\. The feature\-space and second input\-space corrections did not improve validation error beyond the first correction with a single mean and so were not selected; we expect them to contribute more with a larger and more diverse ensemble, and we report the stages that the validation split actually chose\. Two smaller effects are worth isolating\. Reflection test\-time averaging lowers the MLP test error by 0\.11 points, consistent with Proposition[6\.1](https://arxiv.org/html/2609.00389#S6.Thmtheorem1), which guarantees the direction of the change\. Adding the kernel\-flow regularizer \(Section[6\.2](https://arxiv.org/html/2609.00389#S6.SS2)\) to the low\-data MLP did not help here \(it moved test error from 5\.64% to 5\.87%\); Proposition[6\.2](https://arxiv.org/html/2609.00389#S6.Thmtheorem2)says what the term estimates, not that estimating it is useful on every problem, and on this smooth\-output benchmark the plain metric loss was already a good proxy\. Nor does a richer correction kernel help: replacing the single tuned Matérn with a sum of Matérn kernels at several bandwidths leaves the validation error unchanged to two decimals, which is consistent with the residual being close to kernel\-irreducible in the input geometry, and is why the boosted correction stages of Section[3\.3](https://arxiv.org/html/2609.00389#S3.SS3)were not selected\. One positive lesson does emerge from the refiner variants\. Conditioning the refiner on the kernel method’s prediction \(test error 4\.73%\) beats conditioning it on another network’s prediction \(4\.84%\), even though the two conditioning fields have almost the same accuracy; the kernel field carries information complementary to the network’s inductive bias, whereas a second network’s field is largely redundant\. Enlarging the ensemble beyond the residual family is the subject of the next subsection, where the returns of each addition turn out to be predictable in advance from the correlation structure of the members’ errors\.

Table 4:Building the high\-data pipeline\. Mean relativeL2L^\{2\}test error under the 20000\-sample protocol; “\+” rows are cumulative\.StageTest error \(%\)Kernel ridge on loads5\.19Reflection\-averaged neural mean4\.86\+ kernel\-conditioned refiner, stack4\.67\+ diverse members \(Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)\), per\-pixel stack4\.58\+ residual kernel correction4\.55
### 4\.3Architectural diversity and the shared residual

The natural reading of Table[2](https://arxiv.org/html/2609.00389#S4.T2)is that the remaining error is a modeling problem: different architectures make different mistakes, so a more diverse ensemble should stack to a lower number\. We tested this directly\. To the residual family of Section[3](https://arxiv.org/html/2609.00389#S3)we added a Fourier neural operator, a UNet, and a variant of the MLP trained on normalized mean squared error instead of the metric, each trained to the published single\-model level or better \(Table[5](https://arxiv.org/html/2609.00389#S4.T5)\)\. The test refutes the reading\. Figure[2](https://arxiv.org/html/2609.00389#S4.F2)shows the correlation of per\-sample normalized residuals between members: every pair of trained networks, across three architecture families and two losses, correlates between 0\.86 and 0\.96, and the kernel predictor, the most alien member by construction, still correlates 0\.79–0\.89 with all of them\. The members are not making different mistakes\. They are making the same mistake with small private variations\.

The second\-moment identity of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)makes the consequence quantitative before any weights are fitted\. The matrixSSmeasured on the validation split predicts both the optimal convex weights and the achievable stack risk: for the six\-member ensemble the predicted root\-mean\-square relative error of the best mixture is4\.97%4\.97\\%, the fitted stack realizes a mean of4\.61%4\.61\\%, and their ratio is the dispersion factor≈0\.93\\approx 0\.93that is stable across every pipeline we ran\. The weights the prediction assigns fromSSalone, without ever evaluating the metric, concentrated on the refiner and the FNO with a few tenths each on the UNet and the MSE\-trained network, a few percent on the kernel, and zero on the plain MLP whose error is already spanned by the better members of its family, match the weights the direct metric search finds, as part \(i\) of the proposition predicts from the measured\(e1,e2,ϱ\)\(e\_\{1\},e\_\{2\},\\varrho\)\. Stacking does exactly what the correlations permit, and members like these permit little: part \(ii\) puts the infinite\-ensemble floor of an equicorrelated family ate​ϱe\\sqrt\{\\varrho\}, which ate≈4\.7%e\\approx 4\.7\\%and a mean pairwiseϱ≈0\.9\\varrho\\approx 0\.9is about4\.5%4\.5\\%\. Allowing the weights to vary over the grid \(an affine per\-pixel stack, fitted by ridge regression on half the validation split and accepted only because it beat the global weights on the other half\) recovers a little of what global weights cannot see, and the kernel correction adds a few hundredths on top; that is the 4\.55% of Table[2](https://arxiv.org/html/2609.00389#S4.T2), and it is consistent with the floor just computed\.

Where does the common mistake live? Writermr\_\{m\}for membermm’s residual field on a validation case andc=1M​∑mrmc=\\frac\{1\}\{M\}\\sum\_\{m\}r\_\{m\}for the across\-member mean, the component that no amount of uniform averaging removes\. Measured over the validation split and three families \(refiner, FNO, kernel\), the shared component carries90%90\\%of the average residual energy \(Figure[3](https://arxiv.org/html/2609.00389#S4.F3)\)\. Three of its properties matter\. Spatially, its energy concentrates at the two corners of the loaded edge, where the traction boundary condition meets the lateral supports; relative to the local field scale it is three to four times larger there than the grid average, and these are precisely the locations where the finite element solution is least accurate and where interpolation to the41×4141\\times 41grid is most strained\. Spectrally, it is enriched in high radial frequencies relative to the stress fields themselves: the fields put97\.5%97\.5\\%of their energy below radial wavenumber 4, the shared residual only69%69\\%\. And its magnitude is statistically independent of everything we can compute from the input: the correlation of‖c‖\\\|c\\\|with the output norm is0\.010\.01, with the load norm0\.010\.01, with the load’s total variation0\.010\.01\. A component that no architecture avoids, that no input statistic predicts, that lives at the stress concentrations, and that fixed\-scale additive noise would reproduce is most economically read as noise of the data\-generating process itself, finite element and grid\-interpolation error at the singular corners, not as approximation error of the surrogates\.

Two independent measurements corroborate this reading\. First, the residual kernel correction, which regresses the stack’s residual on the load, improves the stack by only a few hundredths of a point atn=19000n=19000whatever the kernel scale: at this sample size the shared residual is not a function of the load in any way a Matérn RKHS onℝ41\\mathbb\{R\}^\{41\}can see\. Second, the error stops responding to data\. Under identical recipes, the MLP’s test error is4\.95%4\.95\\%with35003500training samples,4\.86%4\.86\\%with85008500, and4\.86%4\.86\\%with1900019000\(Figure[5](https://arxiv.org/html/2609.00389#S4.F5)\); the kernel ridge curve still falls, at roughlyn−0\.11n^\{\-0\.11\}, but toward the same region\. A modeling limitation should yield to capacity, diversity, or data\. This error yields to none of them\.

We draw two conclusions\. First, the published plateau of Table[2](https://arxiv.org/html/2609.00389#S4.T2), five methods within two thirds of a point of each other after years of architectural progress, is not evidence that some sixth architecture is missing; every measurement above is what a benchmark whose achievable error is set by its data would produce\. This is an inference from surrogate\-side measurements — certifying it would take regenerated data, mesh\-refinement checks, or repeated solver evaluations, which we did not perform — but on this reading, numbers near4\.5%4\.5\\%sit within a few hundredths of the limit, which is where our pipeline lands \(4\.55%\), and material further progress on this benchmark would come from regenerating the data \(finer meshes near the corners, higher\-order elements, output grids that resolve the concentrations\) rather than from new surrogates\. Second, the practical recipe survives the reinterpretation with its economics clarified: symmetrized accurate means capture what is learnable, stacking buys exactly the decorrelation the members possess \(little, here\), and the kernel correction certifies and mops up the load\-visible remainder\. The components are worth their cost in that order\.

Table 5:Ensemble members, high\-data protocol\. Mean relativeL2L^\{2\}test error with reflection averaging; every member is selected on the validation split\.MemberFamily / lossError \(%\)Residual MLPMLP, metric loss4\.86Residual MLPMLP, normalized MSE4\.71Kernel\-conditioned refinerMLP on\(u,KRR​\(u\)\)\(u,\\mathrm\{KRR\}\(u\)\)4\.73FNOspectral, metric loss4\.70UNetconvolutional, metric loss4\.99Kernel ridge on loadsMatérn\-5/25/25\.19Per\-pixel stack \(validation\-fitted\)–4\.58\+ residual kernel correction–4\.55![Refer to caption](https://arxiv.org/html/2609.00389v1/x2.png)Figure 2:Correlation of per\-sample normalized residuals between ensemble members on the validation split \(test\-set values agree to two decimals\)\. Every pair of trained networks correlates between 0\.86 and 0\.96, the convolutional UNet being the least like the rest; the kernel predictor, the most alien member by construction, still correlates 0\.79–0\.89 with all of them\.![Refer to caption](https://arxiv.org/html/2609.00389v1/x3.png)Figure 3:Anatomy of the shared residualcc\(across\-member mean of the residual fields; refiner, FNO and kernel members\)\. Left to right: rms ofccover the grid; rms of the member\-specific remainders; radial DCT band energies of the stress fields, ofcc, and of the specific parts; and the norm ofccagainst the output norm, whose correlation is0\.010\.01\.
### 4\.4Uncertainty quantification

Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)certifies the corrected surrogate’s error by‖G−m‖K​Pλ​\(u\)\\\|G\-m\\\|\_\{K\}\\,P\_\{\\lambda\}\(u\), and it is tempting to read the design factorPλP\_\{\\lambda\}as a pointwise error bar\. On this benchmark that reading fails, in an instructive way\. Table[6](https://arxiv.org/html/2609.00389#S4.T6)confirms the operator factor behaves as the theory predicts: the computable squared RKHS norm of the fitted correction is about 271×\\timessmaller when the kernel regresses ensemble residuals than when it regresses raw stress fields, so an accurate mean genuinely leaves a smoother remainder and a tighter certificate \(the bound, linear in the norm, by the square root of that\)\. But the design factorPλP\_\{\\lambda\}correlates only weakly with the corrected surrogate’s absolute error \(Pearson 0\.08\), and*negatively*with the benchmark’s relative error\. The cause is measurable:PλP\_\{\\lambda\}is large for loads far from the training set, those loads tend to have large amplitude, large\-amplitude loads produce large\-norm stress fields, and dividing by the output norm makes their relative error small\. IndeedPλP\_\{\\lambda\}correlates 0\.65 with the output norm \(Figure[4](https://arxiv.org/html/2609.00389#S4.F4)\)\. The error is dominated by the neural mean, whose mistakes are not a function of the input\-space geometry thatPλP\_\{\\lambda\}sees, so a purely input\-space power function cannot rank them\.

The obvious alternative fares no better, and with the ensemble of Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)the question can be answered rather than left open\. Ensemble disagreement, the spread of the members’ predictions, is the standard deep\-ensemble error signal; on the full cross\-architecture ensemble its correlation with the corrected surrogate’s absolute error is 0\.08, no better than the 0\.10 of the like\-architecture ensemble\. Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)says why it cannot do better here: disagreement measures the member\-specific components of the error, which carry under a tenth of the residual energy, while the ranking signal for the shared component is invisible to any within\-ensemble statistic, because the members agree precisely where they are jointly wrong\.

What does hold is coverage, for free\. A distribution\-free conformal rescaling, taking the 0\.9\-quantile of‖e‖2/Pλ\\\|e\\\|\_\{2\}/P\_\{\\lambda\}on the validation split and scalingPλP\_\{\\lambda\}by it, yields a band with 91\.6% empirical coverage on the test set at the nominal 90% level\. Conformal validity needs no correlation between the score and the error, only exchangeability, so it survives the confound that defeats the ranking\. One honesty note, expanded in Section[6\.7](https://arxiv.org/html/2609.00389#S6.SS7): the same validation split served model selection, so the exchangeability behind the exact guarantee is approximate here, and the observed coverage should be read as empirical; a calibration split untouched by selection would make the guarantee exact\. The lesson is narrow but real: a valid pointwise bound \(Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)\) is not automatically a useful pointwise indicator, the standard uncertainty signals can both fail on the same problem, and a relative\-error metric has to have its normalization accounted for before a posterior standard deviation means what it appears to\.

Table 6:Squared RKHS norm of the fitted kernel interpolant,tr⁡\(α⊤​K​α\)\\operatorname\{tr\}\(\\alpha^\{\\top\}K\\alpha\), when the same validated kernel and design regress raw stress fields versus ensemble residuals\.![Refer to caption](https://arxiv.org/html/2609.00389v1/x4.png)Figure 4:Left: the Gaussian process posterior standard deviationPλP\_\{\\lambda\}against the corrected surrogate’s relative error on the test set; the association is weak and, because of the output\-amplitude confound, slightly negative\. Right: decile calibration ofPλP\_\{\\lambda\}against absolute error, the theory\-consistent quantity\. A distribution\-free conformal rescaling ofPλP\_\{\\lambda\}attains 91\.6% coverage at the nominal 90% level\.![Refer to caption](https://arxiv.org/html/2609.00389v1/x5.png)Figure 5:Relative test error against training set size for the kernel ridge on loads \(log\-log\), with the full pipeline and the best published results marked in both regimes\.
### 4\.5Spectra and effective dimension

Figure[6](https://arxiv.org/html/2609.00389#S4.F6)shows the spectrum of the Matérn Gram matrix on the loads and the effective dimensiondeff​\(λ\)d\_\{\\mathrm\{eff\}\}\(\\lambda\)computed from it by the exact identity of Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)\. The spectrum decays quickly: a few hundred eigenvalues carry essentially all the mass, the tail falling many orders of magnitude below the leading modes\. At the cross\-validated nugget the effective dimension is a few hundred out ofn=19000n=19000\. This is a statistical statement, not a computational one — the dense factorization costs its fulln3n^\{3\}operations regardless — but it is why the correction generalizes despite interpolating nineteen thousand points: the problem the kernel solves has the statistical dimension of the retained spectrum, not of the sample\. The cross\-validatedλ\\lambdasits at the shoulder of thedeffd\_\{\\mathrm\{eff\}\}curve, past the leading modes and into the rapidly decaying tail, matching the reading of Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)\.

![Refer to caption](https://arxiv.org/html/2609.00389v1/x6.png)Figure 6:Left: eigenvalues of the Matérn Gram matrixK/nK/non the loads \(log scale\); the spectrum decays by many orders of magnitude within a few hundred modes\. Right: the effective dimensiondeff​\(λ\)d\_\{\\mathrm\{eff\}\}\(\\lambda\)from the exact identity of Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11), with the cross\-validated nugget marked\.
### 4\.6Cost and reproducibility

The pipeline was built and run on a single laptop\. The kernel stages are exact: one Cholesky factorization of the19000×1900019000\\times 19000Matérn Gram matrix in double precision takes2424seconds on the laptop’s CPU \(1616threads, OpenBLAS; the matrix itself is2\.92\.9GB\), and prediction is two matrix products; there are no inducing points, stochastic solvers, or preconditioners, and no GPU is needed for the correction\. The factorization performs the fulln3/3≈2\.3×1012n^\{3\}/3\\approx 2\.3\\times 10^\{12\}floating\-point operations — numerical low rank does not reduce the work of a dense Cholesky — and the measured rate, about100100gigaflops per second, is simply what a current multicore CPU delivers; exactness at thisnncosts half a minute, not a cluster\. What the low effective dimension of Section[4\.5](https://arxiv.org/html/2609.00389#S4.SS5)buys is statistical rather than computational: it describes where the solve’s information concentrates at the working nugget, not its flop count\. The neural means are small by the standards of the literature \(a few million parameters\) and train in the metric itself\. All splits are fixed by a published permutation seed, the input is consumed as the 41\-dimensional load after the exact broadcast check of Section[2](https://arxiv.org/html/2609.00389#S2), and every model is selected on the validation split and evaluated once on the test set; we release the splits, the trained means, and the code so the table can be regenerated end to end\.

## 5The OCO\-2 radiative\-transfer emulator

The structural\-mechanics benchmark put the neural mean and the kernel in the balanced regime\. To see the other regime we apply the same components to the radiative\-transfer emulation problem ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\]: per spectral band of the OCO\-2 instrument, learn the map from a reduced atmospheric state \(2020–2424dimensions after the dimension reduction of that paper\) to the reduced radiance \(4040PCA coefficients of the monochromatic spectrum\)\. The OSF project \(u2t8a\) publishes a training pool of2000020000pairs and, separately,20002000test states together with the kernel\-flow emulator’s own predictions on them; we verified the two sets share no point\. We split the pool18000/200018000/2000into training and validation, and every choice in our pipeline — the network checkpoints, the kernel head’s length scale and nugget, and the per\-coordinate winners of the combination below — is made on that validation split; the20002000public test states are used once, for the final numbers\. The comparison is therefore computed on identical test points against the published emulator itself rather than against our reimplementation of it\.

Two error metrics matter, and they disagree in an instructive way\. The*reduced*metric is the relativeL2L^\{2\}error on the4040standardized coefficients\. The*radiance*metric maps predictions back to the monochromatic spectrum through the stored PCA projection and norms before comparing; because the projection is orthogonal, the error*numerator*on the reduced side is exactly the diagonally weighted norm‖sz⊙\(z^−z\)‖\\\|s\_\{z\}\\odot\(\\hat\{z\}\-z\)\\\|, whose weights concentrate almost entirely on the first few coefficients \(the denominator is the full reconstructed\-radiance norm, which adds the reconstruction offsets\)\. It is this diagonal weighting of the residual that the per\-coordinate combination exploits\.

Table 7:OCO\-2 O2 band, scored on the20002000public test states, which are disjoint from the training pool; all model selection happened on a validation split of the pool\. The kernel\-flow row is the emulator ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\], scored from its own stored predictions\. “Features” means the penultimate activations of the trained network; the kernel head is the same exact Matérn solve as everywhere else in this paper\. The two raw\-input rows are the diagnostic\-script fit \(a60006000\-point solve at a single median length scale\), lighter than the feature head’s full\-grid fit; the anisotropy experiment of Section[6\.5](https://arxiv.org/html/2609.00389#S6.SS5), which tunes the raw kernel fully, confirms it stays an order of magnitude above the feature head, so the gap is not an artifact of this\.![Refer to caption](https://arxiv.org/html/2609.00389v1/x7.png)Figure 7:Left: the representation ladder on the O2 band\. The same exact Matérn machinery moves from40%40\\%on the raw input to3\.8%3\.8\\%on the trained network’s features; the kernel\-flow emulator and the network sit in between\. Right: per\-coordinate test RMSE with the reduced coordinates ordered by their radiance weight \(shaded bars\)\. The per\-output\-tuned kernel\-flow emulator is most accurate exactly on the heavily weighted coordinates, the flat\-metric network on the many light ones, and the radiance\-trained network closes the gap where the weight lives; the per\-coordinate combination takes each coordinate from whichever model wins it\.Table[7](https://arxiv.org/html/2609.00389#S5.T7)carries four findings, and Figure[7](https://arxiv.org/html/2609.00389#S5.F7)shows the two central ones\. First, the representation ladder: the same exact kernel machinery scores40%40\\%on the raw input,30%30\\%after rescaling each input coordinate by its measured sensitivity \(the metric term of Proposition[6\.7](https://arxiv.org/html/2609.00389#S6.Thmtheorem7): the effective rank here is1717of2020, so the proposition predicts a constant\-factor gain and no change of rate, the40%→30%40\\%\\\!\\to\\\!30\\%the ARD row records\), and3\.82%3\.82\\%on the trained network’s features, past the network’s own head\. The raw\-input kernel was limited by its features, not by anything kernel\-shaped, and the deep kernel head, an exact solve on learned features, is the strongest single model on the reduced metric — though its margin over the network it feeds on \(3\.82%3\.82\\%against3\.99%3\.99\\%, one seed each\) is small, and we lean on the ladder’s order of magnitude, not on that gap\. Learning the metric rather than setting it by hand does not change this reading\. A per\-dimension metric fit by the kernel\-flows cross\-validation loss ofOwhadi and Yoo \[[17](https://arxiv.org/html/2609.00389#bib.bib7)\]reaches the same30%30\\%on the raw state, and a low\-rank Mahalanobis metric with a thousand more parameters does not improve on it; on the features, where the network has already conditioned the representation, a learned metric gives no gain over the isotropic kernel\. This is what Proposition[6\.7](https://arxiv.org/html/2609.00389#S6.Thmtheorem7)anticipates: the metric is a constant\-factor lever, and the representation is the order\-of\-magnitude one\.

Second, the training metric behaves as a budget\. The flat\-trained and radiance\-trained rows are one architecture and two losses, and each wins the metric it was trained in by a factor of five or more while losing the other\. The mechanism is visible per coordinate: the radiance metric concentrates on the leading coefficient, the kernel\-flow emulator \(whose lengthscales are tuned per output\) predicts that coefficient to0\.00120\.0012root mean square against the flat network’s0\.00540\.0054, while the flat network is three to four times more accurate on the trailing thirty coordinates that the radiance barely sees\. Training the network in the radiance metric closes exactly that gap\. A constant\-elasticity model of this reallocation is too crude — the per\-coordinate exponents scatter widely, and coordinate errors do not decouple because the network shares capacity — but the direction is robust, and the practical rule is plain: train in the metric you report\.

Third, the same\-class floor of the structural\-mechanics study reappears here, one level down\. The residuals of independently seeded copies of the flat network correlate at0\.780\.78on the O2 band and at0\.970\.97and0\.960\.96on WCO2 and SCO2, so part \(ii\) of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)puts the infinite\-seed floors of the flat network at3\.9%3\.9\\%,16\.2%16\.2\\%and8\.0%8\.0\\%, and its seed ensemble sits there\. Passing below the floor took a different kind of member: the Matérn head on learned features reaches3\.82%3\.82\\%,16\.1%16\.1\\%and7\.96%7\.96\\%, and the final combination3\.83%3\.83\\%,16\.1%16\.1\\%and7\.96%7\.96\\%\. Ensembling within an architecture class saturates by the same second\-moment law on both problems; what moved the OCO\-2 numbers an order of magnitude was changing what the members are \(learned features under the kernel\), not how many there are\.

The floor is a property of the architecture class, not of the problem, and it can be probed directly\. Different architectures decorrelate: on WCO2 two residual MLPs of different activation correlate at0\.840\.84, a wide shallow network against the residual MLP at0\.400\.40, and a random\-Fourier\-feature network against it at0\.080\.08to0\.170\.17, all far below the0\.960\.96of two seeds\. Yet the floor did not move, because its other input is accuracy: by part \(i\) of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)a member of errore2e\_\{2\}helps a reference of errore1e\_\{1\}only when their correlation falls belowe1/e2e\_\{1\}/e\_\{2\}, and across a bandwidth sweep the Fourier network was either accurate and correlated \(error31%31\\%, correlation0\.540\.54, just above the0\.530\.53threshold\) or decorrelated and inaccurate \(error above140%140\\%\)\. No architecture we tried was both accurate and decorrelated enough to lower the floor, which is the honest statement of why the residual MLP is the member of record: not that diversity is impossible, but that on this map no more accurate diverse member was found\. Figure[8](https://arxiv.org/html/2609.00389#S5.F8)shows both halves of this: the correlations that make the floor class\-specific, and the accuracy\-decorrelation trade\-off that keeps it in place\. Adding the Fourier members to the per\-coordinate combination changes the reported error by less than a hundredth of a point, in either direction, on all three bands: the honesty\-split blend gives them no weight\.

![Refer to caption](https://arxiv.org/html/2609.00389v1/x8.png)Figure 8:The ensembling floor on WCO2\. Left: residual correlations\. Two seeds of one architecture correlate near the value that sets the floor \(dashed\); different architectures correlate far less, so the floor is a property of the architecture class\. Right: the Fourier\-feature member across bandwidths\. A member helps the ensemble only in the shaded region, below theeref/ϱe\_\{\\mathrm\{ref\}\}/\\varrhocurve of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)\(i\); the accurate settings are too correlated and the decorrelated settings too inaccurate, so none lands there with room to spare\.Fourth, the regimes of Section[6\.4](https://arxiv.org/html/2609.00389#S6.SS4)are decided by these numbers before any stack is fit\. Against the flat network the raw\-input kernel hasϱ=0\.14\\varrho=0\.14anden/ek=0\.10e\_\{n\}/e\_\{k\}=0\.10\(the network’s3\.99%3\.99\\%over the kernel’s40\.47%40\.47\\%\): the condition of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)fails and the kernel is dropped, in contrast to structural mechanics where the same test keeps it\. Because the radiance metric is diagonal in the reduced coordinates, the per\-coordinate combination of Proposition[6\.6](https://arxiv.org/html/2609.00389#S6.Thmtheorem6)is optimal for both metrics simultaneously\. On O2 and SCO2 this yields a single surrogate that beats the kernel\-flow emulator on both metrics at once: O2 reaches3\.83%3\.83\\%against16\.89%16\.89\\%reduced and0\.0267%0\.0267\\%against0\.0448%0\.0448\\%radiance, and SCO2 reaches7\.96%7\.96\\%against16\.14%16\.14\\%reduced and0\.043%0\.043\\%against0\.115%0\.115\\%radiance\. WCO2 is a partial exception, stated plainly\. We did not produce a single combined WCO2 surrogate, and its two best heads split the metrics: the flat\-feature head reaches16\.1%16\.1\\%reduced, below the emulator’s24\.1%24\.1\\%, but0\.101%0\.101\\%radiance,*above*the emulator’s0\.060%0\.060\\%; the radiance\-trained head reaches0\.035%0\.035\\%radiance but at48%48\\%reduced\. On WCO2, then, each metric is won by a different head, and unlike the other two bands no single model beats the emulator on both\. The code releases the per\-band numbers alongside this paper\.

### 5\.1The error is limited by data, and yields to more of it

The structural\-mechanics error was flat in the sample size because it is set by the noise of the finite element data; the emulator error is not, and the contrast is the sharpest practical difference between the two problems\. The left panel of Figure[9](https://arxiv.org/html/2609.00389#S5.F9)fixes the O2 task and varies only the number of training pairs, scoring each on a disjoint held\-out set\. The corrected surrogate’s reduced\-radiance error falls as a clean power law innn, fitted slope−0\.68\-0\.68, from74%74\\%at a few hundred pairs to4\.3%4\.3\\%atn=17700n=17700, still descending at the largest sample size we can afford and approaching the3\.8%3\.8\\%our full pipeline reaches on all1800018000pairs \(Table[7](https://arxiv.org/html/2609.00389#S5.T7)\)\. The exact kernel on the raw state improves far more slowly, so at these sizes the network is the component that turns data into accuracy\. The residual on this task is sample\-limited rather than noise\-limited: unlike the structural\-mechanics floor it recedes as pairs are added, so the largest lever on this problem is simply more of them, and pairs are cheap because the forward model is a program\. The complementary case is the full\-resolution58→304858\\to 3048map, where the error is instead flat innnacross the few hundred retrievals we hold: that task is limited by its representation, not its sample count, and the reduction the emulator applies is what moves it\.

The same reading holds well beyond OCO\-2, on a benchmark large enough to watch the two regimes trade places\. ClimSim\[[24](https://arxiv.org/html/2609.00389#bib.bib33)\]is a climate emulation dataset of about ten million samples that maps a124124\-dimensional atmospheric state to128128physics tendencies, with published neural baselines\. The right panel of Figure[9](https://arxiv.org/html/2609.00389#S5.F9)sweeps the training size across more than three orders of magnitude at a fixed test set, in the coefficient of determination the benchmark reports \(over the outputs with non\-negligible variance\)\. The exact kernel on the state is data\-efficient: itsR2R^\{2\}is already about0\.40\.4at ten thousand samples and holds there\. The neural mean is data\-hungry: badly underfit at a thousand samples, withR2R^\{2\}well below zero, it crosses the kernel near a few hundred thousand samples and keeps rising, reachingR2=0\.56R^\{2\}=0\.56at two million and0\.570\.57at three and a half million, closing on the published multilayer\-perceptron baseline near0\.60\.6\. Which member is the stronger one is therefore a matter of the sample size alone—the kernel below the crossover, the network above it—and the two regimes of Section[6\.4](https://arxiv.org/html/2609.00389#S6.SS4)gain a data axis: the accurate mean is worth building precisely once there is enough data to train it past the kernel it would otherwise defer to\.

![Refer to caption](https://arxiv.org/html/2609.00389v1/x9.png)Figure 9:Test error against training\-set size, at fixed test set\. Left: the OCO\-2 O2 task, reduced\-radiance error of the corrected surrogate; it falls as a power law \(fitted slope−0\.68\-0\.68\) to4\.3%4\.3\\%atn=17700n=17700, approaching the3\.8%3\.8\\%the full pipeline reaches on all1800018000pairs \(dotted\)\. Right: ClimSim, testR2R^\{2\}over the outputs with non\-negligible variance; the exact kernel is flat and data\-efficient while the neural mean, underfit at smallnn, overtakes it near a few hundred thousand samples and closes on the published baseline \(dotted\) by a few million\.

## 6Theory

Throughout this section𝒰=ℝp\\mathcal\{U\}=\\mathbb\{R\}^\{p\}denotes the \(discretized\) input space and𝒱=ℝq\\mathcal\{V\}=\\mathbb\{R\}^\{q\}the output space,G:𝒰→𝒱G:\\mathcal\{U\}\\to\\mathcal\{V\}the target operator,μ\\muthe input distribution, and\(u1,v1\),…,\(uN,vN\)\(u\_\{1\},v\_\{1\}\),\\dots,\(u\_\{N\},v\_\{N\}\)i\.i\.d\. draws withvn=G​\(un\)v\_\{n\}=G\(u\_\{n\}\); the data carry no sampling noise in the usual sense, coming from a deterministic solver, though Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)has something to say about solver error\. We writeℓ​\(w,v\)=‖w−v‖2/‖v‖2\\ell\(w,v\)=\\\|w\-v\\\|\_\{2\}/\\\|v\\\|\_\{2\}for the relative error andR​\(G^\)=𝔼u∼μ​ℓ​\(G^​\(u\),G​\(u\)\)R\(\\widehat\{G\}\)=\\mathbb\{E\}\_\{u\\sim\\mu\}\\,\\ell\(\\widehat\{G\}\(u\),G\(u\)\)for the risk, which is the quantity reported in Section[4](https://arxiv.org/html/2609.00389#S4)\. The results below are largely elementary, but each one licenses a specific design decision in the pipeline of Section[3](https://arxiv.org/html/2609.00389#S3), and each has a measurable footprint in the experiments: one statement per stage, from the symmetry handling and the kernel\-flow term through stacking, the residual correction, its coverage guarantee, and the choice of nugget\. Proofs are collected in Appendix[A](https://arxiv.org/html/2609.00389#A1)\.

### 6\.1Symmetry averaging does not increase the risk

The elasticity problem behind the benchmark is invariant under the reflectionx1↦1−x1x\_\{1\}\\mapsto 1\-x\_\{1\}: the domain, the boundary partition, and the constitutive model are all symmetric, the input distribution has a reflection\-invariant covariance, and reflecting the load reflects the stress field\. On the discrete grid this is expressed by two permutation matrices:S∈ℝp×pS\\in\\mathbb\{R\}^\{p\\times p\}reverses the load samples andT∈ℝq×qT\\in\\mathbb\{R\}^\{q\\times q\}reverses the rows of the stress field\. Both are orthogonal involutions\. Section[2](https://arxiv.org/html/2609.00389#S2)describes a direct data\-driven check of the equivarianceG​\(S​u\)=T​G​\(u\)G\(Su\)=T\\,G\(u\); we treat it as an assumption here\.

###### Proposition 6\.1\(symmetrization\)\.

AssumeG​\(S​u\)=T​G​\(u\)G\(Su\)=TG\(u\)forμ\\mu\-a\.e\.uu, thatμ\\muis invariant underSS, and thatTTis an orthogonal involution \(T2=IT^\{2\}=I\)\. For a measurableG^:𝒰→𝒱\\widehat\{G\}:\\mathcal\{U\}\\to\\mathcal\{V\}define its symmetrization

G^S​\(u\)=12​\(G^​\(u\)\+T​G^​\(S​u\)\)\.\\widehat\{G\}\_\{S\}\(u\)\\;=\\;\\tfrac\{1\}\{2\}\\left\(\\widehat\{G\}\(u\)\+T\\,\\widehat\{G\}\(Su\)\\right\)\.Then, for every lossℓ​\(⋅,v\)\\ell\(\\cdot,v\)that is convex in its first argument and satisfiesℓ​\(T​w,T​v\)=ℓ​\(w,v\)\\ell\(Tw,Tv\)=\\ell\(w,v\),

𝔼u∼μ​ℓ​\(G^S​\(u\),G​\(u\)\)≤𝔼u∼μ​ℓ​\(G^​\(u\),G​\(u\)\)\.\\mathbb\{E\}\_\{u\\sim\\mu\}\\,\\ell\\big\(\\widehat\{G\}\_\{S\}\(u\),G\(u\)\\big\)\\;\\leq\\;\\mathbb\{E\}\_\{u\\sim\\mu\}\\,\\ell\\big\(\\widehat\{G\}\(u\),G\(u\)\\big\)\.

The relative error satisfies both hypotheses \(TTorthogonal givesℓ​\(T​w,T​v\)=ℓ​\(w,v\)\\ell\(Tw,Tv\)=\\ell\(w,v\)\)\. Proposition[6\.1](https://arxiv.org/html/2609.00389#S6.Thmtheorem1)justifies the two uses of the symmetry in the pipeline: averaging each trained model with its reflected evaluation at test time can only help on average, and training on reflection\-augmented data is ordinary empirical risk minimization for the symmetrized class\. The measured effect is small but strictly positive for every model we trained \(Table[4](https://arxiv.org/html/2609.00389#S4.T4)\)\.

### 6\.2What the kernel\-flow regularizer estimates

The kernel\-flow \(KF\) loss of Owhadi and Yoo\[[17](https://arxiv.org/html/2609.00389#bib.bib7)\], in theℓ2\\ell^\{2\}variant used byYoo and Owhadi \[[23](https://arxiv.org/html/2609.00389#bib.bib8)\]to regularize deep networks, enters our ablation study as an optional term in the training of the neural means\. Its motivation is often stated informally \(“a kernel is good if half the data predicts the rest”\)\. The following observation makes the statement exact\. Letzi=ϕθ​\(ui\)z\_\{i\}=\\phi\_\{\\theta\}\(u\_\{i\}\)be the features of a batchB=\{1,…,b\}B=\\\{1,\\dots,b\\\}\(bbeven\), letkkbe a positive definite kernel on feature space, and forC⊂BC\\subset B,\|C\|=b/2\|C\|=b/2, letICI\_\{C\}denote thekk\-interpolant of\{\(zj,vj\):j∈C\}\\\{\(z\_\{j\},v\_\{j\}\):j\\in C\\\}\. The batch KF loss is

e2​\(θ;C\)=∑i∈B‖vi−IC​\(zi\)‖22\.e\_\{2\}\(\\theta;C\)\\;=\\;\\sum\_\{i\\in B\}\\big\\\|v\_\{i\}\-I\_\{C\}\(z\_\{i\}\)\\big\\\|\_\{2\}^\{2\}\.\(1\)
###### Proposition 6\.2\(KF loss is a leave\-half\-out estimate\)\.

LetCCbe drawn uniformly among the subsets ofBBof sizeb/2b/2and assumekkrestricted to\{zi\}i∈B\\\{z\_\{i\}\\\}\_\{i\\in B\}is strictly positive definite\. Thene2​\(θ;C\)=∑i∈B∖C‖vi−IC​\(zi\)‖22e\_\{2\}\(\\theta;C\)=\\sum\_\{i\\in B\\setminus C\}\\\|v\_\{i\}\-I\_\{C\}\(z\_\{i\}\)\\\|\_\{2\}^\{2\}, and

𝔼C​\[e2​\(θ;C\)\]=b2​𝔼C,i∼Unif​\(B∖C\)​\[‖vi−IC​\(zi\)‖22\],\\mathbb\{E\}\_\{C\}\\big\[e\_\{2\}\(\\theta;C\)\\big\]\\;=\\;\\frac\{b\}\{2\}\\;\\mathbb\{E\}\_\{C,\\,i\\,\\sim\\,\\mathrm\{Unif\}\(B\\setminus C\)\}\\big\[\\big\\\|v\_\{i\}\-I\_\{C\}\(z\_\{i\}\)\\big\\\|\_\{2\}^\{2\}\\big\],i\.e\.2b​e2\\tfrac\{2\}\{b\}\\,e\_\{2\}is an unbiased estimator, over the split randomness, of the expected squared error committed on a held\-out point by the kernel predictor trained on a random half of the batch\. When\(B,C\)\(B,C\)are resampled at every step, the stochastic gradient∇θe2​\(θ;C\)\\nabla\_\{\\theta\}e\_\{2\}\(\\theta;C\)is unbiased for∇θ𝔼B,C​\[e2​\(θ;C\)\]\\nabla\_\{\\theta\}\\,\\mathbb\{E\}\_\{B,C\}\[e\_\{2\}\(\\theta;C\)\]whenever the expectation and derivative commute \(e\.g\. under local Lipschitz domination\)\.

The training loss we minimize for a neural mean is a sum of an empirical risk term and, optionally,β​e2\\beta\\,e\_\{2\}: the first term is a linear functional of the empirical measure of the batch and measures fit, while Proposition[6\.2](https://arxiv.org/html/2609.00389#S6.Thmtheorem2)shows the second measures the*generalization of a kernel method built on the learned features*\. This is the sense in which a KF\-regularized network is trained to be a good feature map for the kernel correction that follows\.

### 6\.3Stacking on a held\-out split is safe

Between the means and the correction sits one more fitted object: the convex weights of the stacked ensemble, chosen to minimize the empirical risk on the validation split\. Two elementary facts justify this step\. First, by convexity, a convex combination of predictors is never worse than the corresponding weighted average of their risks, so stacking cannot be hurt by a weak member beyond its weight\. Second, the weights live on a low\-dimensional simplex and are fitted on hundreds of samples, so the selection cannot meaningfully overfit; the following proposition quantifies this with explicit constants\.

###### Proposition 6\.3\(validation stacking\)\.

Letf1,…,fM:𝒰→𝒱f\_\{1\},\\dots,f\_\{M\}:\\mathcal\{U\}\\to\\mathcal\{V\}be fixed maps, letΔ=\{w∈ℝM:wm≥0,∑mwm=1\}\\Delta=\\\{w\\in\\mathbb\{R\}^\{M\}:w\_\{m\}\\geq 0,\\sum\_\{m\}w\_\{m\}=1\\\}, and forw∈Δw\\in\\DeltawriteFw=∑mwm​fmF\_\{w\}=\\sum\_\{m\}w\_\{m\}f\_\{m\}\. Assume the members’ per\-sample relative errors are bounded:ℓ​\(fm​\(u\),G​\(u\)\)≤B\\ell\(f\_\{m\}\(u\),G\(u\)\)\\leq Bforμ\\mu\-a\.e\.uuand everymm\.

1. \(i\)For everyw∈Δw\\in\\Delta,ℓ​\(Fw​\(u\),G​\(u\)\)≤∑mwm​ℓ​\(fm​\(u\),G​\(u\)\)\\ell\(F\_\{w\}\(u\),G\(u\)\)\\leq\\sum\_\{m\}w\_\{m\}\\,\\ell\(f\_\{m\}\(u\),G\(u\)\)pointwise; in particularR​\(Fw\)≤∑mwm​R​\(fm\)≤maxm⁡R​\(fm\)R\(F\_\{w\}\)\\leq\\sum\_\{m\}w\_\{m\}R\(f\_\{m\}\)\\leq\\max\_\{m\}R\(f\_\{m\}\), andℓ​\(Fw​\(u\),G​\(u\)\)≤B\\ell\(F\_\{w\}\(u\),G\(u\)\)\\leq B\.
2. \(ii\)Letw^\\widehat\{w\}minimize the empirical riskR^​\(w\)=1m​∑i=1mℓ​\(Fw​\(u~i\),G​\(u~i\)\)\\widehat\{R\}\(w\)=\\tfrac\{1\}\{m\}\\sum\_\{i=1\}^\{m\}\\ell\(F\_\{w\}\(\\tilde\{u\}\_\{i\}\),G\(\\tilde\{u\}\_\{i\}\)\)overΔ\\Delta, computed on a validation sampleu~1,…,u~m∼μ\\tilde\{u\}\_\{1\},\\dots,\\tilde\{u\}\_\{m\}\\sim\\muindependent off1,…,fMf\_\{1\},\\dots,f\_\{M\}\. Then for everyδ∈\(0,1\)\\delta\\in\(0,1\), with probability at least1−δ1\-\\delta, R​\(Fw^\)≤minw∈Δ⁡R​\(Fw\)\+8​Bm\+B​2​\(\(M−1\)​log⁡\(m​\(M−1\)\+1\)\+log⁡\(2/δ\)\)m\.R\(F\_\{\\widehat\{w\}\}\)\\;\\leq\\;\\min\_\{w\\in\\Delta\}R\(F\_\{w\}\)\\;\+\\;\\frac\{8B\}\{m\}\\;\+\\;B\\sqrt\{\\frac\{2\\big\(\(M\-1\)\\log\\\!\\big\(m\(M\-1\)\+1\\big\)\+\\log\(2/\\delta\)\\big\)\}\{m\}\}\.

WithM=4M=4members,m=1000m=1000validation samples andδ=0\.05\\delta=0\.05the right\-hand excess evaluates to about0\.24​B0\.24\\,B; atB=0\.25B=0\.25, the largest per\-sample member error we observe, this caps the selection cost at roughly six percentage points\. The bound is conservative, as such bounds are: the realized validation\-to\-test gap of the stack in Section[4](https://arxiv.org/html/2609.00389#S4)is two orders of magnitude smaller\. Its role is the scalingM​log⁡m/m\\sqrt\{M\\log m/m\}, which says that fitting a handful of convex weights on a thousand held\-out samples is statistically almost free; this is why the pipeline spends its validation data there and on the kernel hyperparameters and nowhere else\.

### 6\.4What the stack can achieve

Proposition[6\.3](https://arxiv.org/html/2609.00389#S6.Thmtheorem3)says fitting the weights is safe; it does not say mixing is worth anything\. That is decided by a single measurable matrix\. For a mapffwrite the*normalized residual*atuuasρf​\(u\)=\(f​\(u\)−G​\(u\)\)/‖G​\(u\)‖2\\rho\_\{f\}\(u\)=\\big\(f\(u\)\-G\(u\)\\big\)/\\\|G\(u\)\\\|\_\{2\}, so thatℓ​\(f​\(u\),G​\(u\)\)=‖ρf​\(u\)‖2\\ell\(f\(u\),G\(u\)\)=\\\|\\rho\_\{f\}\(u\)\\\|\_\{2\}\.

###### Proposition 6\.4\(second\-moment identity for convex stacks\)\.

Letf1,…,fMf\_\{1\},\\dots,f\_\{M\}be fixed maps and defineS∈ℝM×MS\\in\\mathbb\{R\}^\{M\\times M\}bySm​k=𝔼u∼μ​⟨ρfm​\(u\),ρfk​\(u\)⟩S\_\{mk\}=\\mathbb\{E\}\_\{u\\sim\\mu\}\\,\\langle\\rho\_\{f\_\{m\}\}\(u\),\\rho\_\{f\_\{k\}\}\(u\)\\rangle\. Forwwin the simplexΔM−1\\Delta^\{M\-1\}letfw=∑mwm​fmf\_\{w\}=\\sum\_\{m\}w\_\{m\}f\_\{m\}\. Then

𝔼u∼μ​ℓ​\(fw​\(u\),G​\(u\)\)2=w⊤​S​w,\\mathbb\{E\}\_\{u\\sim\\mu\}\\,\\ell\\big\(f\_\{w\}\(u\),G\(u\)\\big\)^\{2\}\\;=\\;w^\{\\top\}S\\,w,so the squared\-metric risk of every convex stack is determined bySS, and the best achievable isminw∈ΔM−1⁡w⊤​S​w\\min\_\{w\\in\\Delta^\{M\-1\}\}w^\{\\top\}Sw\. In particular:

1. \(i\)\(two members\) ifSShas diagonal\(e12,e22\)\(e\_\{1\}^\{2\},e\_\{2\}^\{2\}\)withe1≤e2e\_\{1\}\\leq e\_\{2\}and correlationϱ=S12/\(e1​e2\)\\varrho=S\_\{12\}/\(e\_\{1\}e\_\{2\}\), mixing in the weaker member strictly helps if and only ifϱ<e1/e2\\varrho<e\_\{1\}/e\_\{2\}, in which case the optimum is interior with valuee12​e22​\(1−ϱ2\)/\(e12\+e22−2​ϱ​e1​e2\)e\_\{1\}^\{2\}e\_\{2\}^\{2\}\(1\-\\varrho^\{2\}\)/\(e\_\{1\}^\{2\}\+e\_\{2\}^\{2\}\-2\\varrho e\_\{1\}e\_\{2\}\);
2. \(ii\)\(equicorrelated family\) ifSm​m=e2S\_\{mm\}=e^\{2\}andSm​k=ϱ​e2S\_\{mk\}=\\varrho\\,e^\{2\}for allm≠km\\neq kwithϱ≥0\\varrho\\geq 0, thenminw⁡w⊤​S​w=e2​\(ϱ\+\(1−ϱ\)/M\)\\min\_\{w\}w^\{\\top\}Sw=e^\{2\}\\big\(\\varrho\+\(1\-\\varrho\)/M\\big\), attained at uniform weights and decreasing toe2​ϱe^\{2\}\\varrhoasM→∞M\\to\\infty\.

Exact equicorrelation never holds in practice, but for convex weights the floor survives one\-sided bounds on the entries ofSS, which are what one actually measures\.

###### Corollary 6\.5\(ensembling floor\)\.

If every member satisfiesSm​m≥e¯2S\_\{mm\}\\geq\\bar\{e\}^\{\\,2\}and every pair satisfiesSm​k≥ϱ¯​e¯2S\_\{mk\}\\geq\\bar\{\\varrho\}\\,\\bar\{e\}^\{\\,2\}withϱ¯∈\[0,1\]\\bar\{\\varrho\}\\in\[0,1\], then every convex combinationfwf\_\{w\}obeys

𝔼​ℓ​\(fw​\(u\),G​\(u\)\)2=w⊤​S​w≥e¯2​\(ϱ¯\+\(1−ϱ¯\)​‖w‖22\)≥e¯2​ϱ¯\.\\mathbb\{E\}\\,\\ell\\big\(f\_\{w\}\(u\),G\(u\)\\big\)^\{2\}\\;=\\;w^\{\\top\}Sw\\;\\geq\\;\\bar\{e\}^\{\\,2\}\\big\(\\bar\{\\varrho\}\+\(1\-\\bar\{\\varrho\}\)\\,\\\|w\\\|\_\{2\}^\{2\}\\big\)\\;\\geq\\;\\bar\{e\}^\{\\,2\}\\,\\bar\{\\varrho\}\.No number of additional members with the same error level and mutual correlations moves the ensemble’s root\-mean\-square relative error belowe¯​ϱ¯\\bar\{e\}\\sqrt\{\\bar\{\\varrho\}\}\.

A second consequence of the same decomposition explains the final OCO\-2 surrogate\. Evaluation metrics there are diagonal in the output coordinates \(the radiance metric is the weightingszs\_\{z\}\), and diagonal metrics decouple\.

###### Proposition 6\.6\(per\-coordinate combination under diagonal metrics\)\.

For weightsv=\(vm​j\)v=\(v\_\{mj\}\)that may depend on the output coordinate, letfv​\(u\)j=∑mvm​j​fm​\(u\)jf\_\{v\}\(u\)\_\{j\}=\\sum\_\{m\}v\_\{mj\}f\_\{m\}\(u\)\_\{j\}\. For any diagonal metric with positive weightsww,

𝔼​‖diag​\(w\)​\(fv​\(u\)−G​\(u\)\)‖22=∑jwj2​𝔼​\(fv​\(u\)j−G​\(u\)j\)2,\\mathbb\{E\}\\big\\\|\\mathrm\{diag\}\(w\)\\big\(f\_\{v\}\(u\)\-G\(u\)\\big\)\\big\\\|\_\{2\}^\{2\}=\\sum\_\{j\}w\_\{j\}^\{2\}\\;\\mathbb\{E\}\\big\(f\_\{v\}\(u\)\_\{j\}\-G\(u\)\_\{j\}\\big\)^\{2\},so the minimizing weights solveqqseparate problems, one per output coordinate, none of which involvesww\. The per\-coordinate optimum is the same for every diagonal metric, and it weakly dominates every combination whose weights are shared across coordinates, for all such metrics simultaneously\.

The practical content: when a problem is scored in several diagonal metrics at once, the combination stage does not need to know which one matters\. Fit the best combination coordinate by coordinate on validation and the result serves all of them; only the members need metric\-aware training, which is where the weighted mean of Section[5](https://arxiv.org/html/2609.00389#S5)enters\.

The corollary is the quantity this paper keeps measuring, at two different levels\. Across architecture families on the structural\-mechanics benchmark the pairwise correlations exceed0\.860\.86at member errors near4\.7%4\.7\\%, a floor of about4\.4%4\.4\\%, and the corrected stack lands within a tenth of it\. Across random seeds of one architecture on the OCO\-2 bands the correlations are0\.780\.78,0\.970\.97and0\.960\.96at member errors of4\.4%4\.4\\%,16\.4%16\.4\\%and8\.2%8\.2\\%, floors of3\.9%3\.9\\%,16\.2%16\.2\\%and8\.0%8\.0\\%, and the measured combinations sit on all three\. When an ensemble is at its floor the corollary says where improvement cannot come from; on OCO\-2 it came from changing the members \(the kernel head on learned features\) rather than adding more of them\. The reported metric is the mean rather than the root mean square of‖ρf‖2\\\|\\rho\_\{f\}\\\|\_\{2\}, and the two differ by the dispersion of the per\-sample error; across every pipeline in this paper their ratio stays within a few percent of0\.930\.93, so predictions made throughSStransfer to the reported metric essentially unchanged\. The proposition frames the diversity experiment of Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3): what a new member is worth is legible in the row it adds toSS, before any weights are fitted\. On this benchmark the rows all look alike, pairwise correlations of 0\.86–0\.96 among trained networks of three architecture families and two losses, 0\.79–0\.89 against the kernel predictor, so the floor sits just below the best single member, and the measured stack lands on it\. Part \(i\) is the same computation specialized to two members; it correctly predicts, from measured\(e1,e2,ϱ\)\(e\_\{1\},e\_\{2\},\\varrho\)alone, which members receive weight in the fitted stack and which are dropped\. Section[5](https://arxiv.org/html/2609.00389#S5)applies the same test on a problem where the answer comes out the other way: there the kernel’s error is inflated by a factor the next subsection quantifies, the thresholde1/e2e\_\{1\}/e\_\{2\}collapses, and the kernel is dropped even at a residual correlation of0\.140\.14\.

### 6\.5Metric and representation: when a fixed kernel loses to a network

The structural\-mechanics benchmark has the neural mean and the isotropic kernel within half a point of each other\. On the OCO\-2 emulation problem of Section[5](https://arxiv.org/html/2609.00389#S5)the same fixed kernel is an order of magnitude behind the same network\. Two mechanisms can produce such a gap, and they separate cleanly\.

Suppose the operator factors asG​\(u\)=g​\(A​u\)G\(u\)=g\(Au\)withA∈ℝd×dA\\in\\mathbb\{R\}^\{d\\times d\}nonsingular andggof unitHsH^\{s\}norm, and thatμ\\muhas a density bounded above and below on a compact domain\. Order the singular valuesσ1≥⋯≥σd\\sigma\_\{1\}\\geq\\dots\\geq\\sigma\_\{d\}ofAA; they are the rates at whichGGvaries along the input directions\. SupposeAAhas*approximate rank*rrat the sample size, meaning the trailing singular values fall below the achievable resolution,σr\+1≲n−1/r\\sigma\_\{r\+1\}\\lesssim n^\{\-1/r\}, so the imageA​𝒰A\\,\\mathcal\{U\}is effectivelyrr\-dimensional at scalen−1/rn^\{\-1/r\}\.

###### Proposition 6\.7\(isotropic and adapted kernel rates\)\.

LetG^niso\\widehat\{G\}^\{\\mathrm\{iso\}\}\_\{n\}be kernel ridge regression with a Matérn kernel of smoothnessssand a validation\-optimal single length scale, andG^nard\\widehat\{G\}^\{\\mathrm\{ard\}\}\_\{n\}the same construction in the metricu↦A​uu\\mapsto Au\. Up to constants depending on\(σi\)\(\\sigma\_\{i\}\),gg,ssand the density bounds, and up to logarithmic factors innn, with high probability overnni\.i\.d\. samples,

𝔼​‖G−G^niso‖≲σ1s​n−s/d,𝔼​‖G−G^nard‖≲n−s/r,\\mathbb\{E\}\\big\\\|G\-\\widehat\{G\}^\{\\mathrm\{iso\}\}\_\{n\}\\big\\\|\\;\\lesssim\\;\\sigma\_\{1\}^\{\\,s\}\\,n^\{\-s/d\},\\qquad\\mathbb\{E\}\\big\\\|G\-\\widehat\{G\}^\{\\mathrm\{ard\}\}\_\{n\}\\big\\\|\\;\\lesssim\\;n^\{\-s/r\},so the ratio of the two bounds is of orderσ1s​ns​\(1/r−1/d\)\\sigma\_\{1\}^\{\\,s\}\\,n^\{\\,s\(1/r\-1/d\)\}\. WhenAAis close to full rank \(r≈dr\\approx d\) the two rates coincide and only the constantσ1s\\sigma\_\{1\}^\{\\,s\}separates them\.

The proof \(Appendix[A](https://arxiv.org/html/2609.00389#A1)\) is the scattered\-data estimate‖G−G^n‖≲hXs​‖G‖Hs\\\|G\-\\widehat\{G\}\_\{n\}\\\|\\lesssim h\_\{X\}^\{\\,s\}\\\|G\\\|\_\{H^\{s\}\}applied to the two designs: one length scale must fill alldddirections, sohX≍n−1/dh\_\{X\}\\asymp n^\{\-1/d\}and the norm carriesσ1s\\sigma\_\{1\}^\{\\,s\}; the adapted metric, whenAAhas a spectral gap, confines the design torrresolved directions and improves the fill distance\. Both displays are upper bounds under the stated idealization — the exact metricAA, a fixedgg, and the length\-scale selection folded into the logarithmic factors — and we claim no matching lower bound; the proposition is used here to calibrate what a metric can and cannot buy, not as a sharp rate theorem\. The proposition bounds the*metric*part of a network’s advantage, the part an anisotropic kernel would recover, and it says this part is a rate gain only under an approximate\-rank gap; without one it is the constantσ1s\\sigma\_\{1\}^\{\\,s\}alone\. The rest is*representational*: a stationary kernel of fixed smoothness reaches only its native space, at any metric, while the network adapts its features to the map\. The two parts are separated experimentally by handing the kernel the anisotropic metric and seeing how much of the gap closes\. On the OCO\-2 band of Section[5](https://arxiv.org/html/2609.00389#S5)the input has effective rank1717of2020, close to full: the proposition predicts the metric can supply only a constant, not a rate, and the measurement agrees,1\.41\.4of a tenfold gap\. The remaining factor is representational, and its direct evidence is that the same kernel machinery applied to the network’s learned features recovers the network’s accuracy and slightly exceeds it\.

The representational term has a precise reading through the same optimal\-recovery bound that Section[6\.6](https://arxiv.org/html/2609.00389#S6.SS6)makes exact\. That bound factors the error as‖G−m‖ℋ⋅Pλ\\\|G\-m\\\|\_\{\\mathcal\{H\}\}\\cdot P\_\{\\lambda\}, a product of the target’s norm in the kernel’s native space and a design factor set by the Gram spectrum\. Changing the features changes both in principle, but on the O2 band only one moves\. The effective dimension of the Matérn Gram, the quantity Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)controls, is nearly identical on the raw input and on the learned features \(atn=4000n=4000and a matched median length scale,39953995against39943994at a10−810^\{\-8\}nugget, with equal leading eigenvalue mass\), so the design factor is essentially fixed\. The native\-space norm of the target is not: interpolating the same outputs through the two kernels, the raw\-input kernel carries the target at squared norm1\.64×1061\.64\\times 10^\{6\}and the feature kernel at3\.87×1043\.87\\times 10^\{4\}, a factor of4242\. The network does not condition the kernel better; it moves the target to a low\-norm corner of a native space of the same size, where the exact solve reaches it\.

This is not an accident of the O2 band; it is what a pulled\-back kernel does\. Writingφ\\varphifor the feature map andkkfor the base kernel, the deep kernel head useskφ​\(x,x′\)=k​\(φ​\(x\),φ​\(x′\)\)k\_\{\\varphi\}\(x,x^\{\\prime\}\)=k\(\\varphi\(x\),\\varphi\(x^\{\\prime\}\)\)\.

###### Proposition 6\.8\(feature pullback\)\.

The native space ofkφk\_\{\\varphi\}isℋkφ=\{f∘φ:f∈ℋk\}\\mathcal\{H\}\_\{k\_\{\\varphi\}\}=\\\{f\\circ\\varphi:f\\in\\mathcal\{H\}\_\{k\}\\\}with‖g‖ℋkφ=min⁡\{‖f‖ℋk:f∘φ=g​on​𝒳\}\\\|g\\\|\_\{\\mathcal\{H\}\_\{k\_\{\\varphi\}\}\}=\\min\\\{\\\|f\\\|\_\{\\mathcal\{H\}\_\{k\}\}:f\\circ\\varphi=g\\text\{ on \}\\mathcal\{X\}\\\}\. In particular, if the target factors through the features asG=h∘φG=h\\circ\\varphiwithh∈ℋkh\\in\\mathcal\{H\}\_\{k\}, then‖G‖ℋkφ≤‖h‖ℋk\\\|G\\\|\_\{\\mathcal\{H\}\_\{k\_\{\\varphi\}\}\}\\leq\\\|h\\\|\_\{\\mathcal\{H\}\_\{k\}\}, and the optimal\-recovery bound of Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)for the feature\-kernel regression is governed by‖h‖ℋk\\\|h\\\|\_\{\\mathcal\{H\}\_\{k\}\}rather than by the norm ofGGin the raw\-input space\.

The proof \(Appendix[A](https://arxiv.org/html/2609.00389#A1)\) is the standard pullback identity for reproducing kernels\. Its content here is the reading: a fixed kernel on the raw input must carry the whole warped mapGG, whose native\-space norm is large whenGGhas fine structure; the same kernel on the features carries only the unwarpedhh, and training drivesφ\\varphitoward exactly the factorization that makeshhsimple\. The measured factor of4242is that norm gap, and it is the representational advantage that an anisotropic metric, which rescales the input but cannot re\-express the map, leaves untouched \(recovering only the factor1\.41\.4\)\.

### 6\.6An error bound for the residual kernel correction

The final stage of the pipeline is kernel ridge regression of the ensemble residual\. The bound below is a vector\-valued, residual form of the classical power\-function estimate for kernel interpolation\[[21](https://arxiv.org/html/2609.00389#bib.bib12)\]and of the optimal\-recovery viewpoint ofOwhadi and Scovel \[[16](https://arxiv.org/html/2609.00389#bib.bib9)\]; we include the short argument in Appendix[A](https://arxiv.org/html/2609.00389#A1)to keep the two factors explicit\. Fix the meanm:𝒰→𝒱m:\\mathcal\{U\}\\to\\mathcal\{V\}\(in the pipeline, the stacked ensemble; in the statement belowmmis any fixed map independent of the correction sample\) and writer=G−mr=G\-mfor the residual operator with componentsrjr\_\{j\},j=1,…,qj=1,\\dots,q\. Letkkbe a positive definite kernel on𝒰\\mathcal\{U\}with RKHSℋk\\mathcal\{H\}\_\{k\}, and assumerj∈ℋkr\_\{j\}\\in\\mathcal\{H\}\_\{k\}for alljjwith

‖r‖K2:=∑j=1q‖rj‖ℋk2<∞\.\\\|r\\\|\_\{K\}^\{2\}:=\\sum\_\{j=1\}^\{q\}\\\|r\_\{j\}\\\|\_\{\\mathcal\{H\}\_\{k\}\}^\{2\}<\\infty\.Given inputsX=\(u1,…,un\)X=\(u\_\{1\},\\dots,u\_\{n\}\), Gram matrixK=k​\(X,X\)K=k\(X,X\)and nuggetλ≥0\\lambda\\geq 0, the correction isr^λ​\(u\)=k​\(u,X\)​\(K\+n​λ​I\)−1​R\\widehat\{r\}\_\{\\lambda\}\(u\)=k\(u,X\)\\,\(K\+n\\lambda I\)^\{\-1\}R,Rn​j=rj​\(un\)R\_\{nj\}=r\_\{j\}\(u\_\{n\}\), and the corrected predictor ism\+r^λm\+\\widehat\{r\}\_\{\\lambda\}\. Let

Pλ​\(u\)2=k​\(u,u\)−k​\(u,X\)​\(K\+n​λ​I\)−1​k​\(X,u\)P\_\{\\lambda\}\(u\)^\{2\}\\;=\\;k\(u,u\)\-k\(u,X\)\\,\(K\+n\\lambda I\)^\{\-1\}k\(X,u\)\(2\)denote the Gaussian process posterior variance with the same nugget \(forλ=0\\lambda=0this is the classical power function of the point setXX\)\.

###### Theorem 6\.9\(certified pointwise bound\)\.

For everyu∈𝒰u\\in\\mathcal\{U\}and everyλ≥0\\lambda\\geq 0,

‖G​\(u\)−m​\(u\)−r^λ​\(u\)‖2≤‖G−m‖K​P~λ​\(u\)≤‖G−m‖K​Pλ​\(u\),\\big\\\|G\(u\)\-m\(u\)\-\\widehat\{r\}\_\{\\lambda\}\(u\)\\big\\\|\_\{2\}\\;\\leq\\;\\\|G\-m\\\|\_\{K\}\\;\\widetilde\{P\}\_\{\\lambda\}\(u\)\\;\\leq\\;\\\|G\-m\\\|\_\{K\}\\;P\_\{\\lambda\}\(u\),whereP~λ​\(u\)2=Pλ​\(u\)2−n​λ​‖\(K\+n​λ​I\)−1​k​\(X,u\)‖22\\widetilde\{P\}\_\{\\lambda\}\(u\)^\{2\}=P\_\{\\lambda\}\(u\)^\{2\}\-n\\lambda\\,\\\|\(K\+n\\lambda I\)^\{\-1\}k\(X,u\)\\\|\_\{2\}^\{2\}\.

Three comments\. First, the inequality is algebraic: it holds pathwise for every fixedmm, every dataset, and everyuu, with no probabilistic assumptions\. In the pipelinemmis itself fit on the same training set; this does not affect the validity of the bound, only the reading of‖G−m‖K\\\|G\-m\\\|\_\{K\}as a random \(realized\) quantity rather than an a priori one\. Second, the bound factorizes into a quantity that depends only on the residual operator \(‖G−m‖K\\\|G\-m\\\|\_\{K\}\) and a quantity that depends only on the kernel and the design \(PλP\_\{\\lambda\}\), and the second factor is exactly the posterior standard deviation that a Gaussian process interpretation of the correction would report\. This is the sense in which the bound is certified: whatever error the corrected surrogate commits atuuis controlled by the reportedPλ​\(u\)P\_\{\\lambda\}\(u\)times a constant that does not depend onuu\. We stress that this is a statement about a valid upper bound, not about the usefulness ofPλP\_\{\\lambda\}as a pointwise error*ranking*\. Section[4\.4](https://arxiv.org/html/2609.00389#S4.SS4)finds that on this benchmarkPλP\_\{\\lambda\}ranks the error poorly, because the neural mean supplies most of the error and the relative metric is confounded by output amplitude; a distribution\-free conformal rescaling ofPλP\_\{\\lambda\}nonetheless attains its nominal coverage on the test set\. Third, the theorem quantifies why one should regress residuals rather than raw targets: the design factorPλP\_\{\\lambda\}is the same in both cases, so the gain of the neural mean is the drop from‖G−v¯‖K\\\|G\-\\bar\{v\}\\\|\_\{K\}\(kernel\-only, meanv¯\\bar\{v\}\) to‖G−m‖K\\\|G\-m\\\|\_\{K\}\. These norms are not directly observable, but the RKHS norms of the*fitted*interpolants,‖r^λ‖K2=tr⁡\(α⊤​K​α\)\\smash\{\\\|\\widehat\{r\}\_\{\\lambda\}\\\|\_\{K\}^\{2\}=\\operatorname\{tr\}\(\\alpha^\{\\top\}K\\alpha\)\}withα=\(K\+n​λ​I\)−1​R\\alpha=\(K\+n\\lambda I\)^\{\-1\}R, are computable lower\-bound proxies, and they drop by a factor of about 271 when the kernel is moved from raw targets to ensemble residuals \(Table[6](https://arxiv.org/html/2609.00389#S4.T6)\)\. The neural mean does not merely reduce the size of the residual inL2L^\{2\}; it leaves behind a residual operator that is genuinely smoother as seen by the kernel\.

### 6\.7Distribution\-free coverage for the reported band

Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)controls the error by‖G−m‖K​Pλ​\(u\)\\\|G\-m\\\|\_\{K\}\\,P\_\{\\lambda\}\(u\), but the constant‖G−m‖K\\\|G\-m\\\|\_\{K\}is not known a priori, and Section[4\.4](https://arxiv.org/html/2609.00389#S4.SS4)shows thatPλP\_\{\\lambda\}alone ranks the error weakly\. What the surrogate reports as uncertainty is therefore notPλP\_\{\\lambda\}itself but a rescaling of it,q​Pλ​\(u\)q\\,P\_\{\\lambda\}\(u\), whose multiplierqqis calibrated on the held\-out split by the split\-conformal rule

q=Q1−α​\(\{si:=‖e​\(u~i\)‖2/Pλ​\(u~i\)\}i=1m\),e​\(u\)=G​\(u\)−m​\(u\)−r^λ​\(u\),q=Q\_\{1\-\\alpha\}\\Big\(\\big\\\{s\_\{i\}:=\\\|e\(\\tilde\{u\}\_\{i\}\)\\\|\_\{2\}/P\_\{\\lambda\}\(\\tilde\{u\}\_\{i\}\)\\big\\\}\_\{i=1\}^\{m\}\\Big\),\\qquad e\(u\)=G\(u\)\-m\(u\)\-\\widehat\{r\}\_\{\\lambda\}\(u\),\(3\)whereQ1−αQ\_\{1\-\\alpha\}is the⌈\(1−α\)​\(m\+1\)⌉\\lceil\(1\-\\alpha\)\(m\+1\)\\rceil\-th smallest value of themmvalidation scores\. The point of this construction is that its coverage needs neither the bound to be tight norPλP\_\{\\lambda\}to be a good error ranking; it needs only exchangeability, which the fixed splits provide\.

###### Proposition 6\.10\(finite\-sample coverage\)\.

Fix the meanmm, the correctionr^λ\\widehat\{r\}\_\{\\lambda\}and the designXX, and letPλ​\(⋅\)\>0P\_\{\\lambda\}\(\\cdot\)\>0\. Suppose the validation and test inputsu~1,…,u~m,u⋆\\tilde\{u\}\_\{1\},\\dots,\\tilde\{u\}\_\{m\},u^\{\\star\}are exchangeable \(in particular i\.i\.d\. fromμ\\mu\) and drawn independently ofm,r^λ,Xm,\\widehat\{r\}\_\{\\lambda\},X\. Letqqbe the conformal multiplier \([3](https://arxiv.org/html/2609.00389#S6.E3)\) and define the bandC​\(u\)=\{v:‖v−m​\(u\)−r^λ​\(u\)‖2≤q​Pλ​\(u\)\}C\(u\)=\\\{v:\\\|v\-m\(u\)\-\\widehat\{r\}\_\{\\lambda\}\(u\)\\\|\_\{2\}\\leq q\\,P\_\{\\lambda\}\(u\)\\\}\. Then

1−α≤ℙ​\(G​\(u⋆\)∈C​\(u⋆\)\)≤1−α\+1m\+1,1\-\\alpha\\;\\leq\\;\\mathbb\{P\}\\big\(G\(u^\{\\star\}\)\\in C\(u^\{\\star\}\)\\big\)\\;\\leq\\;1\-\\alpha\+\\frac\{1\}\{m\+1\},the upper bound holding when the scores are almost surely distinct\.

The two\-sided statement is the standard guarantee of split conformal prediction\[[20](https://arxiv.org/html/2609.00389#bib.bib35),[10](https://arxiv.org/html/2609.00389#bib.bib36)\], transported to the functional\-output setting by taking the nonconformity score to be the relative residual norms=‖e​\(u\)‖2/Pλ​\(u\)s=\\\|e\(u\)\\\|\_\{2\}/P\_\{\\lambda\}\(u\); the proof, a rank argument on the exchangeable scores, is in Appendix[A](https://arxiv.org/html/2609.00389#A1)\. Three points make it the right closing statement for the uncertainty analysis\. It is close to the quantity we compute: the reported coverage of 91\.6% at the nominal1−α=90%1\-\\alpha=90\\%in Section[4\.4](https://arxiv.org/html/2609.00389#S4.SS4)realizes this proposition withm=1000m=1000andα=0\.1\\alpha=0\.1\. The1\.61\.6\-point excess over nominal is not the1/\(m\+1\)1/\(m\+1\)slack, which is only a tenth of a point; it is the finite\-sample fluctuation of coverage on a single test set, whose standard deviationα​\(1−α\)/m≈0\.95\\sqrt\{\\alpha\(1\-\\alpha\)/m\}\\approx 0\.95point places 91\.6% about1\.71\.7standard deviations above nominal, well inside the guarantee\. One caveat on the hypotheses: our calibration scores are computed on the same held\-out split used to select the neural\-mean checkpoints and tune the correction, so the exchangeability the proposition assumes is only approximate here; the selection is low\-complexity and the band over\-covers, but a calibration split disjoint from model selection would make the guarantee exact\. It is agnostic to everything the earlier results leave uncertain: the neural mean may be biased,PλP\_\{\\lambda\}may correlate with the error weakly or with the wrong sign, and the bound of Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)may be loose, yet the band still covers at the stated rate\. And it is where the uncertainty analysis ends: among the candidate uncertainty signals we examine, the conformally rescaledPλP\_\{\\lambda\}is the only one that carries a guarantee, so it is the one the surrogate reports\.

### 6\.8Effective dimension from the Gram spectrum

The remaining design choice is the nuggetλ\\lambda, selected by cross\-validation in the pipeline\. The following exact identity connects that choice to the spectrum of the Gram matrix and, through it, to random matrix descriptions of the design\.

###### Lemma 6\.11\(effective dimension and the empirical Stieltjes transform\)\.

LetKKbe a symmetric positive semidefiniten×nn\\times nmatrix with eigenvaluesλ1,…,λn≥0\\lambda\_\{1\},\\dots,\\lambda\_\{n\}\\geq 0, letm^​\(z\)=1n​∑i=1n\(λi/n−z\)−1\\widehat\{m\}\(z\)=\\tfrac\{1\}\{n\}\\sum\_\{i=1\}^\{n\}\(\\lambda\_\{i\}/n\-z\)^\{\-1\}be the Stieltjes transform of the empirical spectral distribution ofK/nK/n\. Then for everyλ\>0\\lambda\>0the effective dimension of ridge regression with nuggetλ\\lambdasatisfies

deff​\(λ\):=tr⁡\(K​\(K\+n​λ​I\)−1\)=n​\(1−λ​m^​\(−λ\)\)\.d\_\{\\mathrm\{eff\}\}\(\\lambda\)\\;:=\\;\\operatorname\{tr\}\\\!\\big\(K\(K\+n\\lambda I\)^\{\-1\}\\big\)\\;=\\;n\\big\(1\-\\lambda\\,\\widehat\{m\}\(\-\\lambda\)\\big\)\.In particular, if the empirical spectral distribution ofK/nK/nconverges weakly to a limit lawFFasn→∞n\\to\\infty, thendeff​\(λ\)/n→1−λ​mF​\(−λ\)d\_\{\\mathrm\{eff\}\}\(\\lambda\)/n\\to 1\-\\lambda\\,m\_\{F\}\(\-\\lambda\)withmFm\_\{F\}the Stieltjes transform ofFF\.

The identity is unconditional; the random\-matrix content enters through the choice ofFF\. For the benchmark’s simplest reference model, isotropic linear features, the limit can be evaluated in closed form\.

###### Corollary 6\.12\(effective dimension under a Marčenko–Pastur limit\)\.

LetK=X​X⊤K=XX^\{\\top\}withX∈ℝn×pX\\in\\mathbb\{R\}^\{n\\times p\}having i\.i\.d\. entries of mean zero and varianceσ2\\sigma^\{2\}, and letp/n→γ∈\(0,1\]p/n\\to\\gamma\\in\(0,1\]\. Then

deff​\(λ\)n⟶γ​\(1−λ​m​\(−λ\)\),m​\(−λ\)=\(λ\+σ2​\(1−γ\)\)2\+4​γ​σ2​λ−\(λ\+σ2​\(1−γ\)\)2​γ​σ2​λ,\\frac\{d\_\{\\mathrm\{eff\}\}\(\\lambda\)\}\{n\}\\;\\longrightarrow\\;\\gamma\\big\(1\-\\lambda\\,m\(\-\\lambda\)\\big\),\\qquad m\(\-\\lambda\)=\\frac\{\\sqrt\{\\big\(\\lambda\+\\sigma^\{2\}\(1\-\\gamma\)\\big\)^\{2\}\+4\\gamma\\sigma^\{2\}\\lambda\}\-\\big\(\\lambda\+\\sigma^\{2\}\(1\-\\gamma\)\\big\)\}\{2\\gamma\\sigma^\{2\}\\lambda\},almost surely, wheremmis the Stieltjes transform of the Marčenko–Pastur law with ratioγ\\gammaand scaleσ2\\sigma^\{2\}\. In a simulation withn=6000n=6000,p=1800p=1800,σ2=1\.7\\sigma^\{2\}=1\.7, the formula agrees with the empiricaldeff​\(λ\)/nd\_\{\\mathrm\{eff\}\}\(\\lambda\)/nto four decimal places acrossλ∈\[10−3,10\]\\lambda\\in\[10^\{\-3\},10\]\.

Two uses\. Quantitatively, the formula makes the nugget–capacity trade\-off explicit for the bulk of a feature Gram matrix: eigenvalue mass of Marčenko–Pastur type is absorbed or discarded byλ\\lambdaaccording to a single closed\-form curve, and outlying \(spiked\) eigenvaluesxsx\_\{s\}simply add their termsxs/\(xs\+n​λ\)x\_\{s\}/\(x\_\{s\}\+n\\lambda\)on top\. Qualitatively, it separates the two regimes visible in Figure[6](https://arxiv.org/html/2609.00389#S4.F6): the Matérn Gram matrix on the loads is nothing like a Marčenko–Pastur bulk \(its spectrum decays exponentially, which is whydeffd\_\{\\mathrm\{eff\}\}is a few hundred out of 19000\), whereas the Gram matrices of learned penultimate features do show a bulk\-plus\-spikes shape, and for them the corollary describes how much of the bulk the cross\-validated nugget retains\. For the feature Gram matrices of the trained networks we observe the familiar picture of a bulk that is well fitted by a Marčenko–Pastur law\[[14](https://arxiv.org/html/2609.00389#bib.bib14)\]together with a small number of outlying eigenvalues carrying the regression signal, as in spiked covariance models\[[1](https://arxiv.org/html/2609.00389#bib.bib15)\]; for the Matérn Gram matrix on the 41\-dimensional loads the spectrum decays rapidly anddeffd\_\{\\mathrm\{eff\}\}is small \(≈\\approxa few hundred at the cross\-validatedλ\\lambda, out ofn=19000n=19000; Figure[6](https://arxiv.org/html/2609.00389#S4.F6)\)\. This gives a post\-hoc reading of the cross\-validated nugget:λ\\lambdalands wheredeff​\(λ\)d\_\{\\mathrm\{eff\}\}\(\\lambda\)has absorbed the outlying eigenvalues and the leading bulk, and further decreasingλ\\lambdabuys capacity precisely where the spectrum carries little signal\. It also explains why the exact solve is affordable: the correction is numerically a problem of dimensiondeffd\_\{\\mathrm\{eff\}\}, notnn\.

The effective dimension is also exactly the total budget of the design factor from Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)\. WritingA=K\+n​λ​IA=K\+n\\lambda I, the posterior variances at the training inputs sum to

∑i=1nPλ​\(ui\)2=tr⁡\(K\)−tr⁡\(K​A−1​K\)=tr⁡\(n​λ​K​A−1\)=n​λ​deff​\(λ\),\\sum\_\{i=1\}^\{n\}P\_\{\\lambda\}\(u\_\{i\}\)^\{2\}=\\operatorname\{tr\}\(K\)\-\\operatorname\{tr\}\(KA^\{\-1\}K\)=\\operatorname\{tr\}\\\!\\big\(n\\lambda\\,KA^\{\-1\}\\big\)=n\\lambda\\,d\_\{\\mathrm\{eff\}\}\(\\lambda\),usingK​A−1​K=K−n​λ​K​A−1KA^\{\-1\}K=K\-n\\lambda\\,KA^\{\-1\}\. So the samedeffd\_\{\\mathrm\{eff\}\}that reads the nugget off the spectrum also fixes, up to the factorn​λn\\lambda, the total posterior\-variance mass the correction distributes over the design; the two halves of this section are one quantity seen from two sides\.

## 7Discussion

#### What moved the number, and what stopped it\.

The improvement over the published band decomposes unevenly\. Accurate neural means trained on the reported metric and averaged over the reflection symmetry do most of the work; stacking members from different architecture families adds exactly what their measured residual correlations permit, about a tenth of a point here; the kernel corrections contribute a few hundredths more and, with them, the certificate of Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)and the calibrated band of Proposition[6\.10](https://arxiv.org/html/2609.00389#S6.Thmtheorem10)\. We did not find evidence that any published architecture was mis\-tuned by its authors; the components compose, and the composition had not been tried on this benchmark\. What stops further movement is not a missing architecture\. Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)locates the remaining error in a component that all families share, that concentrates where the finite element data are least reliable, that no input statistic predicts, and that neither capacity, nor diversity, nor more data reduces\. The margin over the published methods other than PARA\-Net is real, and the match with PARA\-Net is exact to within run\-to\-run noise, but both should be read for what they are: the last few learnable hundredths above what the evidence of Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)reads as a data\-set floor, not a step on the way to zero\.

#### Relation to neural\-mean Gaussian processes\.

Our correction stage is closest in spirit toMoraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\], who place a neural operator inside the mean of a Gaussian process and fit both jointly\. The differences are operational but they matter at this scale: we fit the mean first and the kernel after \(so the expensive stage is ordinary network training\), we tune the kernel by validation rather than by marginal likelihood, we correct with an exact solve atn=19000n=19000rather than train through the kernel, and we add the feature\-space second stage\. The theory in Section[6](https://arxiv.org/html/2609.00389#S6)applies verbatim to their setting as well; the measured drop in the fitted RKHS norm \(Table[6](https://arxiv.org/html/2609.00389#S4.T6)\) is the quantitative reason residual corrections of accurate means are the right place to spend kernel capacity\.

#### Relation to kernel emulation practice\.

The pipeline is consistent with the experience of the kernel\-flow line of work in emulation at scale\. The closest instance is the forward\-model emulator ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\]for the OCO\-2 CO2retrievals, in which a Gaussian process with a cross\-validation\-learned kernel replaces the full\-physics radiative transfer code within measurement\-error precision; kernel flows have likewise been used to infer convective\-storm structure from passive microwave observations\[[18](https://arxiv.org/html/2609.00389#bib.bib23)\]\. The lesson of that work matches ours: kernels with data\-adapted hyperparameters are extremely effective once the input is presented in its physical parametrization and the target has been reduced to something smooth\. Here the reduction is performed by the neural ensemble rather than by physics, and the kernel\-flow loss itself admits the leave\-half\-out reading of Proposition[6\.2](https://arxiv.org/html/2609.00389#S6.Thmtheorem2)\. The OCO\-2 setting is also the natural next test bed for the residual coupling studied here\. Its training pairs are simulator\-generated state\-to\-radiance maps of exactly the shape this paper exploits, a low\-dimensional physical state mapped to a smooth high\-dimensional output, the mission’s Level 1 and Level 2 products are publicly distributed through NASA’s Earthdata archive, and the emulation setup ofLamminpääet al\.\[[9](https://arxiv.org/html/2609.00389#bib.bib22)\]specifies the state parametrization and sampling design, so a neural mean with a validated kernel correction and a conformal wrapper can be evaluated there against the pure\-kernel emulator without new data collection\.

#### Limitations\.

The benchmark has a one\-dimensional input function and a fixed geometry; the exact kernel solve exploits both, and problems with high\-dimensional input fields would require the usual approximations \(inducing points, random features\) at some cost to the certificates\. The certified bound of Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)controls the error relative to the RKHS norm of the residual operator, a quantity we can only lower\-bound empirically; the conformal calibration we report is the honest, assumption\-light complement\. Training the means on the metric, the reflection steps, and the stacking are benchmark\-agnostic, but the specific error levels reported here should be read as properties of this dataset and protocol\. Finally, our low\-data protocol reproduces that ofMoraet al\.\[[15](https://arxiv.org/html/2609.00389#bib.bib11)\]up to the unavoidable ambiguity of which 1250 samples are used; we use the first 1250 of the training block and report the same 20000\-sample test set, and we release splits and code so the comparison can be audited\.

#### What the second problem settles\.

The OCO\-2 study keeps every component and inverts the balance, and the pair of problems brackets the design space\. When the network and the kernel tie \(structural mechanics\), the coupling pays: stack them, correct the residual, and the certified and conformal machinery rides along\. When the kernel trails by an order of magnitude \(OCO\-2\), the coupling is not the point; the kernel’s value moves inside the network, as an exact head on its learned features, and the pipeline’s job becomes metric management, training each member in the metric it will be scored in and combining per coordinate\. In both regimes the same two measurements decide everything in advance: the pairwise residual correlations, which set the ensembling floor of Corollary[6\.5](https://arxiv.org/html/2609.00389#S6.Thmtheorem5)at whatever level the members share their errors, and the member error ratio, which the second\-moment identity turns into a keep\-or\-drop verdict for each candidate\. Nothing in that protocol is specific to these two datasets\.

#### Outlook\.

Three directions seem worth pursuing\. First, the benchmark itself: regenerating the dataset with refined meshes at the corner singularities, higher\-order elements, or output grids that resolve the concentrations would move the floor that Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)measures, and would return the problem to being a test of surrogates rather than of its own data\. The diagnostic protocol used there, cross\-family residual correlation, the second\-moment prediction of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4), and the scaling\-in\-nncheck, is cheap to run on any benchmark suspected of the same condition\. Second, the same recipe on the remaining benchmarks ofde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]andBatlleet al\.\[[2](https://arxiv.org/html/2609.00389#bib.bib1)\], where the input functions are genuinely high\-dimensional and the interplay between neural means and kernel corrections should look different\. Third, the uncertainty side:PλP\_\{\\lambda\}is a design quantity, so it can be optimized, and the connection of Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)between the nugget, the spectrum, and the effective dimension suggests principled ways to spend a fixed computational budget on the correction stage\.

## Acknowledgments and funding

The research presented in this paper was supported by the European Research Council \(ERC\) under the European Union’s Horizon 2022 research and innovation programme \(grant agreement No\. 101041711\), by the Simons Foundation as part of the Collaboration on the Mathematical and Scientific Foundations of Deep Learning, by Heights Labs, by the Israel Science Foundation \(grant number 2258/19\), by the Israel Science Foundation \(ISF Grant 4101/25\), and by the U\.S\. National Science Foundation \(NSF Grant OISE\-2401227\)\.

## Declaration of competing interest

The author declares no competing interests\.

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

During the preparation of this work the author made substantial use of generative AI tools, and describes their role here\. A large part of the software was written with Fable 5 \(Anthropic\): the implementation of the experiments and diagnostics, the downloading, assembly and preparation of the benchmark and OCO\-2 datasets, and the adaptation of the publicly released emulator code of Lamminpää et al\. — originally split between Julia and a ReFRACtor\-driving Python layer — into a single Python pipeline for sampling states and reducing radiances\. The main contributions are the author’s: the residual coupling of neural means and exact kernel corrections, the two\-regime framing that organizes the paper, and the supporting theory, all developed with AI assistance in calculation, drafting, and exposition\. Each theorem and its proof was checked independently and more than once — by Fable 5, by ChatGPT Pro \(OpenAI\), and by the author — and the reported numbers were verified against the released per\-run data\. The author reviewed all of the above and takes full responsibility for the content of the manuscript\.

## Data and code availability

The code, the per\-run summaries behind every reported number, and the exact configuration of each experiment are available at[https://github\.com/yspennstate/neural\-means\-kernel\-corrections](https://github.com/yspennstate/neural-means-kernel-corrections)\. The structural\-mechanics, Helmholtz, Navier–Stokes and advection data are distributed through the data record accompanyingde Hoopet al\.\[[5](https://arxiv.org/html/2609.00389#bib.bib2)\]\(data\.caltech\.edu, record20091\); the OCO\-2 emulation data and the kernel\-flow emulator’s predictions through the OSF projectu2t8a; and the ClimSim data through the LEAP repositorysubsampled\_low\_res\.

## References

- \[1\]\(2005\)Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices\.Annals of Probability33\(5\),pp\. 1643–1697\.Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1),[§6\.8](https://arxiv.org/html/2609.00389#S6.SS8.p3.13)\.
- \[2\]P\. Batlle, M\. Darcy, B\. Hosseini, and H\. Owhadi\(2024\)Kernel methods are competitive for operator learning\.Journal of Computational Physics496,pp\. 112549\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§1](https://arxiv.org/html/2609.00389#S1.p2.1),[§2\.1](https://arxiv.org/html/2609.00389#S2.SS1.p1.13),[§2\.2](https://arxiv.org/html/2609.00389#S2.SS2.p1.1),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1),[§4\.1](https://arxiv.org/html/2609.00389#S4.SS1.p1.1),[Table 2](https://arxiv.org/html/2609.00389#S4.T2.5.6.6.1),[Table 3](https://arxiv.org/html/2609.00389#S4.T3.4.3.2.1),[§7](https://arxiv.org/html/2609.00389#S7.SS0.SSS0.Px6.p1.2)\.
- \[3\]A\. Caponnetto and E\. De Vito\(2007\)Optimal rates for the regularized least\-squares algorithm\.Foundations of Computational Mathematics7\(3\),pp\. 331–368\.Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1)\.
- \[4\]M\. Darcy, B\. Hamzi, G\. Livieri, H\. Owhadi, and P\. Tavallali\(2023\)One\-shot learning of stochastic differential equations with data adapted kernels\.Physica D: Nonlinear Phenomena444,pp\. 133583\.Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1)\.
- \[5\]M\. V\. de Hoop, D\. Z\. Huang, E\. Qian, and A\. M\. Stuart\(2022\)The cost\-accuracy trade\-off in operator learning with neural networks\.Journal of Machine Learning1\(3\),pp\. 299–341\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§2\.1](https://arxiv.org/html/2609.00389#S2.SS1.p1.13),[§2\.2](https://arxiv.org/html/2609.00389#S2.SS2.p1.1),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1),[§3\.1](https://arxiv.org/html/2609.00389#S3.SS1.p1.1),[§4\.1](https://arxiv.org/html/2609.00389#S4.SS1.p1.1),[Table 2](https://arxiv.org/html/2609.00389#S4.T2.5.2.2.1),[Table 2](https://arxiv.org/html/2609.00389#S4.T2.5.3.3.1),[Table 2](https://arxiv.org/html/2609.00389#S4.T2.5.4.4.1),[Table 2](https://arxiv.org/html/2609.00389#S4.T2.5.5.5.1),[§7](https://arxiv.org/html/2609.00389#S7.SS0.SSS0.Px6.p1.2),[Data and code availability](https://arxiv.org/html/2609.00389#Sx4.p1.1)\.
- \[6\]A\. Dosovitskiy, L\. Beyer, A\. Kolesnikov, D\. Weissenborn, X\. Zhai, T\. Unterthiner, M\. Dehghani, M\. Minderer, G\. Heigold, S\. Gelly, J\. Uszkoreit, and N\. Houlsby\(2021\)An image is worth 16x16 words: transformers for image recognition at scale\.InInternational Conference on Learning Representations,Cited by:[§3\.1](https://arxiv.org/html/2609.00389#S3.SS1.p3.1)\.
- \[7\]V\. Koltchinskii and E\. Giné\(2000\)Random matrix approximation of spectra of integral operators\.Bernoulli6\(1\),pp\. 113–167\.Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1)\.
- \[8\]N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar\(2023\)Neural operator: learning maps between function spaces with applications to PDEs\.Journal of Machine Learning Research24\(89\),pp\. 1–97\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1)\.
- \[9\]O\. Lamminpää, J\. Susiluoto, J\. Hobbs, J\. McDuffie, A\. Braverman, and H\. Owhadi\(2025\)Forward model emulator for atmospheric radiative transfer using Gaussian processes and cross validation\.Atmospheric Measurement Techniques18,pp\. 673–694\.Cited by:[Table 1](https://arxiv.org/html/2609.00389#S1.T1),[Table 1](https://arxiv.org/html/2609.00389#S1.T1.6.3),[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§1](https://arxiv.org/html/2609.00389#S1.p3.2),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1),[Table 7](https://arxiv.org/html/2609.00389#S5.T7),[Table 7](https://arxiv.org/html/2609.00389#S5.T7.4.2),[Table 7](https://arxiv.org/html/2609.00389#S5.T7.7.4.3.1),[§5](https://arxiv.org/html/2609.00389#S5.p1.7),[§7](https://arxiv.org/html/2609.00389#S7.SS0.SSS0.Px3.p1.1)\.
- \[10\]J\. Lei, M\. G’Sell, A\. Rinaldo, R\. J\. Tibshirani, and L\. Wasserman\(2018\)Distribution\-free predictive inference for regression\.Journal of the American Statistical Association113\(523\),pp\. 1094–1111\.Cited by:[§6\.7](https://arxiv.org/html/2609.00389#S6.SS7.p2.10)\.
- \[11\]Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar\(2021\)Fourier neural operator for parametric partial differential equations\.InInternational Conference on Learning Representations,Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1),[§3\.1](https://arxiv.org/html/2609.00389#S3.SS1.p3.1)\.
- \[12\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis\(2021\)Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.Nature Machine Intelligence3,pp\. 218–229\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1)\.
- \[13\]L\. Lu, X\. Meng, S\. Cai, Z\. Mao, S\. Goswami, Z\. Zhang, and G\. E\. Karniadakis\(2022\)A comprehensive and fair comparison of two neural operators \(with practical extensions\) based on FAIR data\.Computer Methods in Applied Mechanics and Engineering393,pp\. 114778\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1)\.
- \[14\]V\. A\. Marčenko and L\. A\. Pastur\(1967\)Distribution of eigenvalues for some sets of random matrices\.Mathematics of the USSR\-Sbornik1\(4\),pp\. 457–483\.Cited by:[§A\.11](https://arxiv.org/html/2609.00389#A1.SS11.p1.12),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1),[§6\.8](https://arxiv.org/html/2609.00389#S6.SS8.p3.13)\.
- \[15\]C\. Mora, A\. Yousefpour, S\. Hosseinmardi, H\. Owhadi, and R\. Bostanabad\(2025\)Operator learning with Gaussian processes\.Computer Methods in Applied Mechanics and Engineering434,pp\. 117581\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§2\.1](https://arxiv.org/html/2609.00389#S2.SS1.p1.13),[§2\.2](https://arxiv.org/html/2609.00389#S2.SS2.p1.1),[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p1.1),[Figure 1](https://arxiv.org/html/2609.00389#S4.F1),[Figure 1](https://arxiv.org/html/2609.00389#S4.F1.3.2),[Table 3](https://arxiv.org/html/2609.00389#S4.T3),[Table 3](https://arxiv.org/html/2609.00389#S4.T3.3.2),[Table 3](https://arxiv.org/html/2609.00389#S4.T3.4.4.3.1),[Table 3](https://arxiv.org/html/2609.00389#S4.T3.4.6.5.1),[§7](https://arxiv.org/html/2609.00389#S7.SS0.SSS0.Px2.p1.1),[§7](https://arxiv.org/html/2609.00389#S7.SS0.SSS0.Px4.p1.1)\.
- \[16\]H\. Owhadi and C\. Scovel\(2019\)Operator\-adapted wavelets, fast solvers, and numerical homogenization\.Cambridge University Press\.Cited by:[§1](https://arxiv.org/html/2609.00389#S1.p1.1),[§6\.6](https://arxiv.org/html/2609.00389#S6.SS6.p1.10)\.
- \[17\]H\. Owhadi and G\. R\. Yoo\(2019\)Kernel flows: from learning kernels from data into the abyss\.Journal of Computational Physics389,pp\. 22–47\.Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1),[§5](https://arxiv.org/html/2609.00389#S5.p3.9),[§6\.2](https://arxiv.org/html/2609.00389#S6.SS2.p1.10)\.
- \[18\]S\. Prasanth, Z\. S\. Haddad, J\. Susiluoto, A\. J\. Braverman, H\. Owhadi, B\. Hamzi, S\. M\. Hristova\-Veleva, and J\. Turk\(2021\)Kernel flows to infer the structure of convective storms from satellite passive microwave observations\.Note:AGU Fall Meeting Abstracts, abstract A55F\-1445Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1),[§7](https://arxiv.org/html/2609.00389#S7.SS0.SSS0.Px3.p1.1)\.
- \[19\]A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, Ł\. Kaiser, and I\. Polosukhin\(2017\)Attention is all you need\.InAdvances in Neural Information Processing Systems,Vol\.30\.Cited by:[§3\.1](https://arxiv.org/html/2609.00389#S3.SS1.p3.1)\.
- \[20\]V\. Vovk, A\. Gammerman, and G\. Shafer\(2005\)Algorithmic learning in a random world\.Springer\.Cited by:[§6\.7](https://arxiv.org/html/2609.00389#S6.SS7.p2.10)\.
- \[21\]H\. Wendland\(2004\)Scattered data approximation\.Cambridge University Press\.Cited by:[§A\.8](https://arxiv.org/html/2609.00389#A1.SS8.p1.9),[§6\.6](https://arxiv.org/html/2609.00389#S6.SS6.p1.10)\.
- \[22\]D\. H\. Wolpert\(1992\)Stacked generalization\.Neural Networks5\(2\),pp\. 241–259\.Cited by:[§3\.2](https://arxiv.org/html/2609.00389#S3.SS2.p1.4)\.
- \[23\]G\. R\. Yoo and H\. Owhadi\(2021\)Deep regularization and direct training of the inner layers of neural networks with kernel flows\.Physica D: Nonlinear Phenomena426,pp\. 132952\.Cited by:[§2\.4](https://arxiv.org/html/2609.00389#S2.SS4.p2.1),[§6\.2](https://arxiv.org/html/2609.00389#S6.SS2.p1.10)\.
- \[24\]S\. Yu, W\. M\. Hannah, L\. Peng, J\. Lin, M\. A\. Bhouri, R\. Gupta,et al\.\(2023\)ClimSim: a large multi\-scale dataset for hybrid physics\-ML climate emulation\.InAdvances in Neural Information Processing Systems \(Datasets and Benchmarks Track\),Vol\.36\.Cited by:[§5\.1](https://arxiv.org/html/2609.00389#S5.SS1.p2.8)\.

## Appendix AProofs

### A\.1Proof of Proposition[6\.1](https://arxiv.org/html/2609.00389#S6.Thmtheorem1)

By convexity ofℓ​\(⋅,v\)\\ell\(\\cdot,v\),

ℓ​\(G^S​\(u\),G​\(u\)\)≤12​ℓ​\(G^​\(u\),G​\(u\)\)\+12​ℓ​\(T​G^​\(S​u\),G​\(u\)\)\.\\ell\\big\(\\widehat\{G\}\_\{S\}\(u\),G\(u\)\\big\)\\;\\leq\\;\\tfrac\{1\}\{2\}\\,\\ell\\big\(\\widehat\{G\}\(u\),G\(u\)\\big\)\+\\tfrac\{1\}\{2\}\\,\\ell\\big\(T\\widehat\{G\}\(Su\),G\(u\)\\big\)\.SinceTTis an involution, the equivarianceG​\(S​u\)=T​G​\(u\)G\(Su\)=TG\(u\)applied atuugivesT​G​\(S​u\)=T2​G​\(u\)=G​\(u\)TG\(Su\)=T^\{2\}G\(u\)=G\(u\), hence

ℓ​\(T​G^​\(S​u\),G​\(u\)\)=ℓ​\(T​G^​\(S​u\),T​G​\(S​u\)\)=ℓ​\(G^​\(S​u\),G​\(S​u\)\),\\ell\\big\(T\\widehat\{G\}\(Su\),G\(u\)\\big\)=\\ell\\big\(T\\widehat\{G\}\(Su\),TG\(Su\)\\big\)=\\ell\\big\(\\widehat\{G\}\(Su\),G\(Su\)\\big\),using theTT\-invariance of the loss\. Taking expectations and using theSS\-invariance ofμ\\mu\(so thatS​u∼μSu\\sim\\muwhenu∼μu\\sim\\mu\),

𝔼​ℓ​\(G^S​\(u\),G​\(u\)\)≤12​𝔼​ℓ​\(G^​\(u\),G​\(u\)\)\+12​𝔼​ℓ​\(G^​\(S​u\),G​\(S​u\)\)=𝔼​ℓ​\(G^​\(u\),G​\(u\)\)\.∎\\mathbb\{E\}\\,\\ell\\big\(\\widehat\{G\}\_\{S\}\(u\),G\(u\)\\big\)\\leq\\tfrac\{1\}\{2\}\\,\\mathbb\{E\}\\,\\ell\\big\(\\widehat\{G\}\(u\),G\(u\)\\big\)\+\\tfrac\{1\}\{2\}\\,\\mathbb\{E\}\\,\\ell\\big\(\\widehat\{G\}\(Su\),G\(Su\)\\big\)=\\mathbb\{E\}\\,\\ell\\big\(\\widehat\{G\}\(u\),G\(u\)\\big\)\.\\qed

### A\.2Proof of Proposition[6\.2](https://arxiv.org/html/2609.00389#S6.Thmtheorem2)

Strict positive definiteness of the restricted kernel matrix makesICI\_\{C\}the \(unique\) interpolant of the pairs indexed byCC, soIC​\(zj\)=vjI\_\{C\}\(z\_\{j\}\)=v\_\{j\}forj∈Cj\\in Cand the terms of \([1](https://arxiv.org/html/2609.00389#S6.E1)\) indexed byCCvanish, which is the first claim\. For the second, write

e2​\(θ;C\)=∑i∈B∖Cg​\(C,i\),g​\(C,i\):=‖vi−IC​\(zi\)‖22\.e\_\{2\}\(\\theta;C\)\\;=\\;\\sum\_\{i\\in B\\setminus C\}g\(C,i\),\\qquad g\(C,i\):=\\\|v\_\{i\}\-I\_\{C\}\(z\_\{i\}\)\\\|\_\{2\}^\{2\}\.Conditionally onCC, the sum has exactlyb/2b/2terms, soe2​\(θ;C\)=b2​𝔼i∼Unif​\(B∖C\)​\[g​\(C,i\)\]e\_\{2\}\(\\theta;C\)=\\tfrac\{b\}\{2\}\\,\\mathbb\{E\}\_\{i\\sim\\mathrm\{Unif\}\(B\\setminus C\)\}\[g\(C,i\)\], and taking the expectation overCCproves the identity\. The gradient claim is immediate: for fixed\(B,C\)\(B,C\)the mapθ↦e2​\(θ;C\)\\theta\\mapsto e\_\{2\}\(\\theta;C\)is differentiable wherever the Cholesky factorization inICI\_\{C\}is \(the kernel matrix stays positive definite in a neighborhood\), and under a local integrable Lipschitz bound the derivative passes under the expectation, so𝔼B,C​\[∇θe2\]=∇θ𝔼B,C​\[e2\]\\mathbb\{E\}\_\{B,C\}\[\\nabla\_\{\\theta\}e\_\{2\}\]=\\nabla\_\{\\theta\}\\,\\mathbb\{E\}\_\{B,C\}\[e\_\{2\}\]\. ∎

### A\.3Proof of Proposition[6\.3](https://arxiv.org/html/2609.00389#S6.Thmtheorem3)

\(i\) Since∑mwm=1\\sum\_\{m\}w\_\{m\}=1,Fw​\(u\)−G​\(u\)=∑mwm​\(fm​\(u\)−G​\(u\)\)F\_\{w\}\(u\)\-G\(u\)=\\sum\_\{m\}w\_\{m\}\\,\(f\_\{m\}\(u\)\-G\(u\)\), and the triangle inequality gives‖Fw​\(u\)−G​\(u\)‖2≤∑mwm​‖fm​\(u\)−G​\(u\)‖2\\\|F\_\{w\}\(u\)\-G\(u\)\\\|\_\{2\}\\leq\\sum\_\{m\}w\_\{m\}\\\|f\_\{m\}\(u\)\-G\(u\)\\\|\_\{2\}\. Dividing by‖G​\(u\)‖2\\\|G\(u\)\\\|\_\{2\}yields the pointwise claim; taking expectations gives the risk inequality, and bounding eachℓ​\(fm​\(u\),G​\(u\)\)\\ell\(f\_\{m\}\(u\),G\(u\)\)byBBgivesℓ​\(Fw​\(u\),G​\(u\)\)≤B\\ell\(F\_\{w\}\(u\),G\(u\)\)\\leq B\.

\(ii\) Writeℓu​\(w\)=ℓ​\(Fw​\(u\),G​\(u\)\)\\ell\_\{u\}\(w\)=\\ell\(F\_\{w\}\(u\),G\(u\)\)\. First,ℓu\\ell\_\{u\}is Lipschitz onΔ\\Deltafor theℓ1\\ell^\{1\}norm: forw,w′∈Δw,w^\{\\prime\}\\in\\Delta, the reverse triangle inequality and the argument of \(i\) give

\|ℓu​\(w\)−ℓu​\(w′\)\|≤‖∑m\(wm−wm′\)​\(fm​\(u\)−G​\(u\)\)‖2‖G​\(u\)‖2≤B​‖w−w′‖1\.\|\\ell\_\{u\}\(w\)\-\\ell\_\{u\}\(w^\{\\prime\}\)\|\\;\\leq\\;\\frac\{\\big\\\|\\sum\_\{m\}\(w\_\{m\}\-w^\{\\prime\}\_\{m\}\)\(f\_\{m\}\(u\)\-G\(u\)\)\\big\\\|\_\{2\}\}\{\\\|G\(u\)\\\|\_\{2\}\}\\;\\leq\\;B\\,\\\|w\-w^\{\\prime\}\\\|\_\{1\}\.Next, discretize the simplex\. Fork∈ℕk\\in\\mathbb\{N\}let𝒢k=\{w∈Δ:k​w∈ℤM\}\\mathcal\{G\}\_\{k\}=\\\{w\\in\\Delta:\\,kw\\in\\mathbb\{Z\}^\{M\}\\\}; its cardinality is the number of compositions ofkkintoMMnonnegative parts,\(k\+M−1M−1\)≤\(k\+1\)M−1\\binom\{k\+M\-1\}\{M\-1\}\\leq\(k\+1\)^\{M\-1\}\. Givenw∈Δw\\in\\Delta, setwi′=⌊k​wi⌋/kw^\{\\prime\}\_\{i\}=\\lfloor kw\_\{i\}\\rfloor/kfori<Mi<MandwM′=1−∑i<Mwi′w^\{\\prime\}\_\{M\}=1\-\\sum\_\{i<M\}w^\{\\prime\}\_\{i\}; thenw′∈𝒢kw^\{\\prime\}\\in\\mathcal\{G\}\_\{k\}\(the last coordinate is a multiple of1/k1/kand is≥wM≥0\\geq w\_\{M\}\\geq 0because the firstM−1M\-1coordinates only decreased\), and

‖w−w′‖1=∑i<M\(wi−wi′\)\+\(wM′−wM\)=2​∑i<M\(wi−wi′\)≤2​\(M−1\)k\.\\\|w\-w^\{\\prime\}\\\|\_\{1\}=\\sum\_\{i<M\}\(w\_\{i\}\-w^\{\\prime\}\_\{i\}\)\+\\big\(w^\{\\prime\}\_\{M\}\-w\_\{M\}\\big\)=2\\sum\_\{i<M\}\(w\_\{i\}\-w^\{\\prime\}\_\{i\}\)\\;\\leq\\;\\frac\{2\(M\-1\)\}\{k\}\.By \(i\) the per\-sample loss lies in\[0,B\]\[0,B\], so for each fixedw′w^\{\\prime\}Hoeffding’s inequality givesℙ​\(\|R^​\(w′\)−R​\(w′\)\|\>t\)≤2​e−2​m​t2/B2\\mathbb\{P\}\\big\(\|\\widehat\{R\}\(w^\{\\prime\}\)\-R\(w^\{\\prime\}\)\|\>t\\big\)\\leq 2e^\{\-2mt^\{2\}/B^\{2\}\}, and a union bound over𝒢k\\mathcal\{G\}\_\{k\}shows that with probability at least1−δ1\-\\delta,

supw′∈𝒢k\|R^​\(w′\)−R​\(w′\)\|≤B​log⁡\(2​\|𝒢k\|/δ\)2​m\.\\sup\_\{w^\{\\prime\}\\in\\mathcal\{G\}\_\{k\}\}\\big\|\\widehat\{R\}\(w^\{\\prime\}\)\-R\(w^\{\\prime\}\)\\big\|\\;\\leq\\;B\\sqrt\{\\frac\{\\log\\\!\\big\(2\|\\mathcal\{G\}\_\{k\}\|/\\delta\\big\)\}\{2m\}\}\.On this event, for everyw∈Δw\\in\\Delta, approximating by its grid point and using the Lipschitz property for bothRRandR^\\widehat\{R\},

\|R^​\(w\)−R​\(w\)\|≤B​log⁡\(2​\|𝒢k\|/δ\)2​m\+4​B​\(M−1\)k\.\\big\|\\widehat\{R\}\(w\)\-R\(w\)\\big\|\\;\\leq\\;B\\sqrt\{\\frac\{\\log\(2\|\\mathcal\{G\}\_\{k\}\|/\\delta\)\}\{2m\}\}\+\\frac\{4B\(M\-1\)\}\{k\}\.Ifw⋆w^\{\\star\}minimizesRRoverΔ\\Delta, the empirical minimizer satisfiesR​\(Fw^\)≤R^​\(w^\)\+supΔ\|R^−R\|≤R^​\(w⋆\)\+supΔ\|R^−R\|≤R​\(Fw⋆\)\+2​supΔ\|R^−R\|R\(F\_\{\\widehat\{w\}\}\)\\leq\\widehat\{R\}\(\\widehat\{w\}\)\+\\sup\_\{\\Delta\}\|\\widehat\{R\}\-R\|\\leq\\widehat\{R\}\(w^\{\\star\}\)\+\\sup\_\{\\Delta\}\|\\widehat\{R\}\-R\|\\leq R\(F\_\{w^\{\\star\}\}\)\+2\\sup\_\{\\Delta\}\|\\widehat\{R\}\-R\|\. Choosingk=m​\(M−1\)k=m\(M\-1\)and using\|𝒢k\|≤\(m​\(M−1\)\+1\)M−1\|\\mathcal\{G\}\_\{k\}\|\\leq\(m\(M\-1\)\+1\)^\{M\-1\}gives

R​\(Fw^\)≤minw∈Δ⁡R​\(Fw\)\+8​Bm\+B​2​\(\(M−1\)​log⁡\(m​\(M−1\)\+1\)\+log⁡\(2/δ\)\)m,R\(F\_\{\\widehat\{w\}\}\)\\;\\leq\\;\\min\_\{w\\in\\Delta\}R\(F\_\{w\}\)\+\\frac\{8B\}\{m\}\+B\\sqrt\{\\frac\{2\\big\(\(M\-1\)\\log\(m\(M\-1\)\+1\)\+\\log\(2/\\delta\)\\big\)\}\{m\}\},which is the claim\. ∎

### A\.4Proof of Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)

Since∑mwm=1\\sum\_\{m\}w\_\{m\}=1, the stack’s normalized residual isρfw​\(u\)=∑mwm​ρfm​\(u\)\\rho\_\{f\_\{w\}\}\(u\)=\\sum\_\{m\}w\_\{m\}\\rho\_\{f\_\{m\}\}\(u\), hence

𝔼​ℓ​\(fw​\(u\),G​\(u\)\)2=𝔼​‖∑mwm​ρfm​\(u\)‖22=∑m,kwm​wk​𝔼​⟨ρfm​\(u\),ρfk​\(u\)⟩=w⊤​S​w\.\\mathbb\{E\}\\,\\ell\\big\(f\_\{w\}\(u\),G\(u\)\\big\)^\{2\}=\\mathbb\{E\}\\Big\\\|\\sum\_\{m\}w\_\{m\}\\rho\_\{f\_\{m\}\}\(u\)\\Big\\\|\_\{2\}^\{2\}=\\sum\_\{m,k\}w\_\{m\}w\_\{k\}\\,\\mathbb\{E\}\\,\\langle\\rho\_\{f\_\{m\}\}\(u\),\\rho\_\{f\_\{k\}\}\(u\)\\rangle=w^\{\\top\}Sw\.For \(i\), parameterizew=\(1−t,t\)w=\(1\-t,t\);q​\(t\)=w⊤​S​wq\(t\)=w^\{\\top\}Swis a convex quadratic inttwithq′​\(0\)=2​\(ϱ​e1​e2−e12\)q^\{\\prime\}\(0\)=2\(\\varrho e\_\{1\}e\_\{2\}\-e\_\{1\}^\{2\}\), negative precisely whenϱ<e1/e2\\varrho<e\_\{1\}/e\_\{2\}\. In that case the unconstrained minimizer lies in\(0,1\)\(0,1\)and gives the stated interior value; otherwiseq′​\(0\)≥0q^\{\\prime\}\(0\)\\geq 0, the minimum over\[0,1\]\[0,1\]is att=0t=0with valuee12e\_\{1\}^\{2\}, and the second member is dropped\. For \(ii\),S=e2​\(\(1−ϱ\)​I\+ϱ​11⊤\)S=e^\{2\}\\big\(\(1\-\\varrho\)I\+\\varrho\\,\\mathbf\{1\}\\mathbf\{1\}^\{\\top\}\\big\), sow⊤​S​w=e2​\(\(1−ϱ\)​‖w‖22\+ϱ\)w^\{\\top\}Sw=e^\{2\}\\big\(\(1\-\\varrho\)\\\|w\\\|\_\{2\}^\{2\}\+\\varrho\\big\)on the simplex, minimized by the uniform weights, where‖w‖22=1/M\\\|w\\\|\_\{2\}^\{2\}=1/M\. ∎

### A\.5Proof of Corollary[6\.5](https://arxiv.org/html/2609.00389#S6.Thmtheorem5)

Convex weights are nonnegative, so each term ofw⊤​S​w=∑mwm2​Sm​m\+∑m≠kwm​wk​Sm​kw^\{\\top\}Sw=\\sum\_\{m\}w\_\{m\}^\{2\}S\_\{mm\}\+\\sum\_\{m\\neq k\}w\_\{m\}w\_\{k\}S\_\{mk\}can be bounded below entrywise:

w⊤​S​w≥e¯2​∑mwm2\+ϱ¯​e¯2​∑m≠kwm​wk=e¯2​\(‖w‖22\+ϱ¯​\(1−‖w‖22\)\),w^\{\\top\}Sw\\;\\geq\\;\\bar\{e\}^\{\\,2\}\\sum\_\{m\}w\_\{m\}^\{2\}\+\\bar\{\\varrho\}\\,\\bar\{e\}^\{\\,2\}\\sum\_\{m\\neq k\}w\_\{m\}w\_\{k\}=\\bar\{e\}^\{\\,2\}\\Big\(\\\|w\\\|\_\{2\}^\{2\}\+\\bar\{\\varrho\}\\,\(1\-\\\|w\\\|\_\{2\}^\{2\}\)\\Big\),using∑m≠kwm​wk=\(∑mwm\)2−‖w‖22=1−‖w‖22\\sum\_\{m\\neq k\}w\_\{m\}w\_\{k\}=\(\\sum\_\{m\}w\_\{m\}\)^\{2\}\-\\\|w\\\|\_\{2\}^\{2\}=1\-\\\|w\\\|\_\{2\}^\{2\}\. The bracket is a convex combination of11andϱ¯\\bar\{\\varrho\}, hence at leastϱ¯\\bar\{\\varrho\}\. ∎

### A\.6Proof of Proposition[6\.6](https://arxiv.org/html/2609.00389#S6.Thmtheorem6)

The metric is diagonal, so the squared norm splits over coordinates and thejj\-th summandwj2​𝔼​\(∑mvm​j​fm​\(u\)j−G​\(u\)j\)2w\_\{j\}^\{2\}\\,\\mathbb\{E\}\\big\(\\sum\_\{m\}v\_\{mj\}f\_\{m\}\(u\)\_\{j\}\-G\(u\)\_\{j\}\\big\)^\{2\}depends on the weights only through the columnv⋅jv\_\{\\cdot j\}\. Minimizing the sum is thereforeqqindependent minimizations, and the positive factorwj2w\_\{j\}^\{2\}does not move any of theqqargmins, so the optimizer is the same for every positiveww\. A combination with shared weights is the special casevm​j=vmv\_\{mj\}=v\_\{m\}for alljj, a subset of the feasible set of each coordinate problem, which gives the domination\. ∎

### A\.7Proof of Proposition[6\.8](https://arxiv.org/html/2609.00389#S6.Thmtheorem8)

LetT:ℋk→ℝ𝒳T:\\mathcal\{H\}\_\{k\}\\to\\mathbb\{R\}^\{\\mathcal\{X\}\}be the composition mapT​f=f∘φTf=f\\circ\\varphi\. Forx∈𝒳x\\in\\mathcal\{X\}the reproducing property gives\(T​f\)​\(x\)=f​\(φ​\(x\)\)=⟨f,k​\(⋅,φ​\(x\)\)⟩ℋk\(Tf\)\(x\)=f\(\\varphi\(x\)\)=\\langle f,k\(\\cdot,\\varphi\(x\)\)\\rangle\_\{\\mathcal\{H\}\_\{k\}\}, so evaluation ofT​fTfatxxis a bounded functional, and the imageℋkφ:=T​\(ℋk\)\\mathcal\{H\}\_\{k\_\{\\varphi\}\}:=T\(\\mathcal\{H\}\_\{k\}\)carries the quotient norm‖g‖ℋkφ=min⁡\{‖f‖ℋk:T​f=g\}\\\|g\\\|\_\{\\mathcal\{H\}\_\{k\_\{\\varphi\}\}\}=\\min\\\{\\\|f\\\|\_\{\\mathcal\{H\}\_\{k\}\}:Tf=g\\\}, the minimum over the affine subspaceT−1​\(g\)T^\{\-1\}\(g\)\(nonempty exactly whenggis a pullback\)\. This normed space is an RKHS: its reproducing kernel iskφ​\(x,x′\)=⟨k​\(⋅,φ​\(x\)\),k​\(⋅,φ​\(x′\)\)⟩ℋk=k​\(φ​\(x\),φ​\(x′\)\)k\_\{\\varphi\}\(x,x^\{\\prime\}\)=\\langle k\(\\cdot,\\varphi\(x\)\),k\(\\cdot,\\varphi\(x^\{\\prime\}\)\)\\rangle\_\{\\mathcal\{H\}\_\{k\}\}=k\(\\varphi\(x\),\\varphi\(x^\{\\prime\}\)\), since the minimum\-norm representer of evaluation atxxisk​\(⋅,φ​\(x\)\)k\(\\cdot,\\varphi\(x\)\)\. Takingf=hf=hin the minimum gives‖h∘φ‖ℋkφ≤‖h‖ℋk\\\|h\\circ\\varphi\\\|\_\{\\mathcal\{H\}\_\{k\_\{\\varphi\}\}\}\\leq\\\|h\\\|\_\{\\mathcal\{H\}\_\{k\}\}, and substituting this bound into Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)applied with the kernelkφk\_\{\\varphi\}replaces‖G‖\\\|G\\\|by‖h‖ℋk\\\|h\\\|\_\{\\mathcal\{H\}\_\{k\}\}\. ∎

### A\.8Proof of Proposition[6\.7](https://arxiv.org/html/2609.00389#S6.Thmtheorem7)

Both displays are the native\-space estimate for kernel ridge regression: for a Matérn kernel of smoothnesssson a bounded domain, with the nugget chosen as in the pipeline, an estimatorG^n\\widehat\{G\}\_\{n\}of a targetFFin the native space obeys𝔼​‖F−G^n‖≤C​hXs​‖F‖𝒩\\mathbb\{E\}\\\|F\-\\widehat\{G\}\_\{n\}\\\|\\leq C\\,h\_\{X\}^\{\\,s\}\\,\\\|F\\\|\_\{\\mathcal\{N\}\}, withhXh\_\{X\}the fill distance of the design in the kernel’s metric\[[21](https://arxiv.org/html/2609.00389#bib.bib12)\], and fornnsamples from a density bounded above and below on a set of effective dimensiond′d^\{\\prime\}, quasi\-uniformity giveshX≍\(log⁡n/n\)1/d′h\_\{X\}\\asymp\(\\log n/n\)^\{1/d^\{\\prime\}\}with high probability, which we write asn−1/d′n^\{\-1/d^\{\\prime\}\}up to the logarithmic factor absorbed into the statement\.

For the isotropic kernel the metric is Euclidean on the fulldd\-dimensional domain, sohX≍n−1/dh\_\{X\}\\asymp n^\{\-1/d\}; the target isF=g\(A⋅\)F=g\(A\\,\\cdot\), whose Matérn native norm is equivalent to itsHsH^\{s\}norm\. Changing variablesy=A​uy=Augives∥g\(A⋅\)∥Hs≍\|detA\|−1/2σ1s\\\|g\(A\\cdot\)\\\|\_\{H^\{s\}\}\\asymp\|\\det A\|^\{\-1/2\}\\,\\sigma\_\{1\}^\{\\,s\}, since each derivative of order up tossbrings down at most a factorσ1\\sigma\_\{1\}and the Jacobian contributes\|detA\|−1/2\|\\det A\|^\{\-1/2\}inL2L^\{2\}; the\|detA\|−1/2\|\\det A\|^\{\-1/2\}is one of the singular\-value constants absorbed into the statement, leaving the displayedσ1s\\sigma\_\{1\}^\{\\,s\}\. For the adapted kernelkA​\(u,u′\)=k​\(A​u,A​u′\)k\_\{A\}\(u,u^\{\\prime\}\)=k\(Au,Au^\{\\prime\}\)the estimator is the isotropic estimator of the unit\-normggon the design\{A​ui\}\\\{Au\_\{i\}\\\}\. Here the approximate\-rank hypothesis enters: the designA​𝒰A\\,\\mathcal\{U\}has extentσr\+1≲n−1/r\\sigma\_\{r\+1\}\\lesssim n^\{\-1/r\}in each of the trailingd−rd\-rdirections, below the fill distance those directions would otherwise demand, so at resolutionn−1/rn^\{\-1/r\}the design isrr\-dimensional and its fill distance ishX≍n−1/rh\_\{X\}\\asymp n^\{\-1/r\}\(a further\(∏i≤rσi\)\(\\prod\_\{i\\leq r\}\\sigma\_\{i\}\)constant, again absorbed\)\. Substituting the two fill distances into the native\-space estimate and dividing gives the ratio\. If insteadAAis close to full rank,σr\+1\\sigma\_\{r\+1\}is not belown−1/rn^\{\-1/r\}, the trailing directions are resolved, and the adapted fill distance isn−1/dn^\{\-1/d\}like the isotropic one: the rates coincide and onlyσ1s\\sigma\_\{1\}^\{\\,s\}separates the bounds\. ∎

### A\.9Proof of Theorem[6\.9](https://arxiv.org/html/2609.00389#S6.Thmtheorem9)

Fixuuand writeκ=k​\(X,u\)∈ℝn\\kappa=k\(X,u\)\\in\\mathbb\{R\}^\{n\},A=K\+n​λ​IA=K\+n\\lambda I,w=A−1​κw=A^\{\-1\}\\kappa\. For each componentjj, membershiprj∈ℋkr\_\{j\}\\in\\mathcal\{H\}\_\{k\}and the reproducing property give

rj​\(u\)−r^λ,j​\(u\)=⟨rj,k​\(u,⋅\)−∑i=1nwi​k​\(ui,⋅\)⟩ℋk,r\_\{j\}\(u\)\-\\widehat\{r\}\_\{\\lambda,j\}\(u\)=\\Big\\langle r\_\{j\},\\;k\(u,\\cdot\)\-\\sum\_\{i=1\}^\{n\}w\_\{i\}\\,k\(u\_\{i\},\\cdot\)\\Big\\rangle\_\{\\mathcal\{H\}\_\{k\}\},becauser^λ,j​\(u\)=κ⊤​A−1​rj​\(X\)=∑iwi​rj​\(ui\)\\widehat\{r\}\_\{\\lambda,j\}\(u\)=\\kappa^\{\\top\}A^\{\-1\}r\_\{j\}\(X\)=\\sum\_\{i\}w\_\{i\}\\,r\_\{j\}\(u\_\{i\}\)andrj​\(ui\)=⟨rj,k​\(ui,⋅\)⟩r\_\{j\}\(u\_\{i\}\)=\\langle r\_\{j\},k\(u\_\{i\},\\cdot\)\\rangle\. Cauchy–Schwarz yields

\|rj​\(u\)−r^λ,j​\(u\)\|≤‖rj‖ℋk​‖k​\(u,⋅\)−∑iwi​k​\(ui,⋅\)‖ℋk,\\big\|r\_\{j\}\(u\)\-\\widehat\{r\}\_\{\\lambda,j\}\(u\)\\big\|\\;\\leq\\;\\\|r\_\{j\}\\\|\_\{\\mathcal\{H\}\_\{k\}\}\\,\\Big\\\|k\(u,\\cdot\)\-\\sum\_\{i\}w\_\{i\}k\(u\_\{i\},\\cdot\)\\Big\\\|\_\{\\mathcal\{H\}\_\{k\}\},and the second factor is independent ofjj; call itρ​\(u\)\\rho\(u\)\. Expanding,

ρ​\(u\)2=k​\(u,u\)−2​w⊤​κ\+w⊤​K​w=k​\(u,u\)−κ⊤​A−1​\(2​A−K\)​A−1​κ\.\\rho\(u\)^\{2\}=k\(u,u\)\-2\\,w^\{\\top\}\\kappa\+w^\{\\top\}Kw=k\(u,u\)\-\\kappa^\{\\top\}A^\{\-1\}\\big\(2A\-K\\big\)A^\{\-1\}\\kappa\.Since2​A−K=A\+n​λ​I2A\-K=A\+n\\lambda I,

ρ​\(u\)2=k​\(u,u\)−κ⊤​A−1​κ−n​λ​‖A−1​κ‖22=P~λ​\(u\)2≤Pλ​\(u\)2\.\\rho\(u\)^\{2\}=k\(u,u\)\-\\kappa^\{\\top\}A^\{\-1\}\\kappa\-n\\lambda\\,\\\|A^\{\-1\}\\kappa\\\|\_\{2\}^\{2\}=\\widetilde\{P\}\_\{\\lambda\}\(u\)^\{2\}\\;\\leq\\;P\_\{\\lambda\}\(u\)^\{2\}\.Summing the squared componentwise bounds,

‖r​\(u\)−r^λ​\(u\)‖22=∑j\(rj​\(u\)−r^λ,j​\(u\)\)2≤ρ​\(u\)2​∑j‖rj‖ℋk2=P~λ​\(u\)2​‖r‖K2\.∎\\big\\\|r\(u\)\-\\widehat\{r\}\_\{\\lambda\}\(u\)\\big\\\|\_\{2\}^\{2\}=\\sum\_\{j\}\\big\(r\_\{j\}\(u\)\-\\widehat\{r\}\_\{\\lambda,j\}\(u\)\\big\)^\{2\}\\leq\\rho\(u\)^\{2\}\\sum\_\{j\}\\\|r\_\{j\}\\\|^\{2\}\_\{\\mathcal\{H\}\_\{k\}\}=\\widetilde\{P\}\_\{\\lambda\}\(u\)^\{2\}\\,\\\|r\\\|\_\{K\}^\{2\}\.\\qed
Note that the argument does not use interpolation: it holds for everyλ≥0\\lambda\\geq 0, with the \(slightly sharper\) factorP~λ\\widetilde\{P\}\_\{\\lambda\}showing that smoothing can only tighten this particular bound relative to the posterior standard deviationPλP\_\{\\lambda\}that we report\.

### A\.10Proof of Proposition[6\.10](https://arxiv.org/html/2609.00389#S6.Thmtheorem10)

Writesi=‖e​\(u~i\)‖2/Pλ​\(u~i\)s\_\{i\}=\\\|e\(\\tilde\{u\}\_\{i\}\)\\\|\_\{2\}/P\_\{\\lambda\}\(\\tilde\{u\}\_\{i\}\)for the validation scores ands⋆=‖e​\(u⋆\)‖2/Pλ​\(u⋆\)s^\{\\star\}=\\\|e\(u^\{\\star\}\)\\\|\_\{2\}/P\_\{\\lambda\}\(u^\{\\star\}\)for the test score\. Becausemm,r^λ\\widehat\{r\}\_\{\\lambda\},XXandPλP\_\{\\lambda\}are fixed and the inputsu~1,…,u~m,u⋆\\tilde\{u\}\_\{1\},\\dots,\\tilde\{u\}\_\{m\},u^\{\\star\}are exchangeable and independent of them, the scoress1,…,sm,s⋆s\_\{1\},\\dots,s\_\{m\},s^\{\\star\}are exchangeable\. The eventG​\(u⋆\)∈C​\(u⋆\)G\(u^\{\\star\}\)\\in C\(u^\{\\star\}\)is exactly\{s⋆≤q\}\\\{s^\{\\star\}\\leq q\\\}, whereq=s\(⌈\(1−α\)​\(m\+1\)⌉\)q=s\_\{\(\\lceil\(1\-\\alpha\)\(m\+1\)\\rceil\)\}is thek:=⌈\(1−α\)​\(m\+1\)⌉k:=\\lceil\(1\-\\alpha\)\(m\+1\)\\rceil\-th order statistic of the validation scores\.

Consider the augmented samples1,…,sm,s⋆s\_\{1\},\\dots,s\_\{m\},s^\{\\star\}of sizem\+1m\+1and letrank​\(s⋆\)\\mathrm\{rank\}\(s^\{\\star\}\)be its rank \(ties broken uniformly at random, which only helps\)\. By exchangeability the rank is uniform on\{1,…,m\+1\}\\\{1,\\dots,m\+1\\\}, soℙ​\(rank​\(s⋆\)≤k\)=k/\(m\+1\)\\mathbb\{P\}\(\\mathrm\{rank\}\(s^\{\\star\}\)\\leq k\)=k/\(m\+1\)\. Ifs⋆s^\{\\star\}is among thekksmallest of them\+1m\+1values then at mostk−1k\-1of thesis\_\{i\}are below it, sos⋆≤s\(k\)=qs^\{\\star\}\\leq s\_\{\(k\)\}=q; converselys⋆≤qs^\{\\star\}\\leq qforcesrank​\(s⋆\)≤k\\mathrm\{rank\}\(s^\{\\star\}\)\\leq kexcept possibly through ties\. Hence

ℙ​\(s⋆≤q\)≥ℙ​\(rank​\(s⋆\)≤k\)=km\+1=⌈\(1−α\)​\(m\+1\)⌉m\+1≥1−α,\\mathbb\{P\}\(s^\{\\star\}\\leq q\)\\;\\geq\\;\\mathbb\{P\}\(\\mathrm\{rank\}\(s^\{\\star\}\)\\leq k\)=\\frac\{k\}\{m\+1\}=\\frac\{\\lceil\(1\-\\alpha\)\(m\+1\)\\rceil\}\{m\+1\}\\;\\geq\\;1\-\\alpha,which is the lower bound\. When the scores are almost surely distinct there are no ties,\{s⋆≤q\}=\{rank​\(s⋆\)≤k\}\\\{s^\{\\star\}\\leq q\\\}=\\\{\\mathrm\{rank\}\(s^\{\\star\}\)\\leq k\\\}exactly, andk/\(m\+1\)<\(\(1−α\)​\(m\+1\)\+1\)/\(m\+1\)=1−α\+1/\(m\+1\)k/\(m\+1\)<\(\(1\-\\alpha\)\(m\+1\)\+1\)/\(m\+1\)=1\-\\alpha\+1/\(m\+1\), giving the upper bound\. ∎

### A\.11Proof of Corollary[6\.12](https://arxiv.org/html/2609.00389#S6.Thmtheorem12)

The nonzero eigenvalues ofK/n=X​X⊤/nK/n=XX^\{\\top\}\\\!/ncoincide with those of the sample covariance matrixS=X⊤​X/n∈ℝp×pS=X^\{\\top\}X/n\\in\\mathbb\{R\}^\{p\\times p\}, and zero eigenvalues contribute nothing todeff​\(λ\)=∑iλi​\(K\)/\(λi​\(K\)\+n​λ\)d\_\{\\mathrm\{eff\}\}\(\\lambda\)=\\sum\_\{i\}\\lambda\_\{i\}\(K\)/\(\\lambda\_\{i\}\(K\)\+n\\lambda\), so

deff​\(λ\)n=pn⋅1p​∑j=1pλj​\(S\)λj​\(S\)\+λ=pn​\(1−λ​m^S​\(−λ\)\),\\frac\{d\_\{\\mathrm\{eff\}\}\(\\lambda\)\}\{n\}=\\frac\{p\}\{n\}\\cdot\\frac\{1\}\{p\}\\sum\_\{j=1\}^\{p\}\\frac\{\\lambda\_\{j\}\(S\)\}\{\\lambda\_\{j\}\(S\)\+\\lambda\}=\\frac\{p\}\{n\}\\Big\(1\-\\lambda\\,\\widehat\{m\}\_\{S\}\(\-\\lambda\)\\Big\),by the computation of Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)applied toSS, withm^S\\widehat\{m\}\_\{S\}the Stieltjes transform of the empirical spectral distribution ofSS\. By the Marčenko–Pastur theorem\[[14](https://arxiv.org/html/2609.00389#bib.bib14)\], this distribution converges weakly, almost surely, to the law with ratioγ\\gammaand scaleσ2\\sigma^\{2\}, and sincex↦\(x\+λ\)−1x\\mapsto\(x\+\\lambda\)^\{\-1\}is bounded and continuous on\[0,∞\)\[0,\\infty\),m^S​\(−λ\)→m​\(−λ\)\\widehat\{m\}\_\{S\}\(\-\\lambda\)\\to m\(\-\\lambda\)\. It remains to evaluatem​\(−λ\)m\(\-\\lambda\)\. The Stieltjes transform of the Marčenko–Pastur law satisfies the self\-consistent equation

m​\(z\)=1σ2​\(1−γ−γ​z​m​\(z\)\)−z,m\(z\)=\\frac\{1\}\{\\sigma^\{2\}\\big\(1\-\\gamma\-\\gamma z\\,m\(z\)\\big\)\-z\},i\.e\.γ​σ2​z​m2\+\(z−σ2​\(1−γ\)\)​m\+1=0\\gamma\\sigma^\{2\}z\\,m^\{2\}\+\\big\(z\-\\sigma^\{2\}\(1\-\\gamma\)\\big\)m\+1=0\. Atz=−λz=\-\\lambdathe discriminant is\(λ\+σ2​\(1−γ\)\)2\+4​γ​σ2​λ\>0\\big\(\\lambda\+\\sigma^\{2\}\(1\-\\gamma\)\\big\)^\{2\}\+4\\gamma\\sigma^\{2\}\\lambda\>0, and of the two real roots

m​\(−λ\)=σ2​\(1−γ\)\+λ±\(λ\+σ2​\(1−γ\)\)2\+4​γ​σ2​λ−2​γ​σ2​λm\(\-\\lambda\)=\\frac\{\\sigma^\{2\}\(1\-\\gamma\)\+\\lambda\\pm\\sqrt\{\\big\(\\lambda\+\\sigma^\{2\}\(1\-\\gamma\)\\big\)^\{2\}\+4\\gamma\\sigma^\{2\}\\lambda\}\}\{\-2\\gamma\\sigma^\{2\}\\lambda\}only the one with the minus sign in the numerator is positive, asm​\(−λ\)=∫\(x\+λ\)−1​𝑑F​\(x\)\>0m\(\-\\lambda\)=\\int\(x\+\\lambda\)^\{\-1\}dF\(x\)\>0requires; rearranging gives the stated expression\. ∎

### A\.12Proof of Lemma[6\.11](https://arxiv.org/html/2609.00389#S6.Thmtheorem11)

DiagonalizeKKwith eigenvaluesλi\\lambda\_\{i\}\. Then

tr⁡\(K​\(K\+n​λ​I\)−1\)=∑i=1nλiλi\+n​λ=∑i=1n\(1−n​λλi\+n​λ\)=n−λ​∑i=1n1λi/n\+λ=n​\(1−λ​m^​\(−λ\)\),\\operatorname\{tr\}\\big\(K\(K\+n\\lambda I\)^\{\-1\}\\big\)=\\sum\_\{i=1\}^\{n\}\\frac\{\\lambda\_\{i\}\}\{\\lambda\_\{i\}\+n\\lambda\}=\\sum\_\{i=1\}^\{n\}\\Big\(1\-\\frac\{n\\lambda\}\{\\lambda\_\{i\}\+n\\lambda\}\\Big\)=n\-\\lambda\\sum\_\{i=1\}^\{n\}\\frac\{1\}\{\\lambda\_\{i\}/n\+\\lambda\}=n\\big\(1\-\\lambda\\,\\widehat\{m\}\(\-\\lambda\)\\big\),usingm^​\(−λ\)=1n​∑i\(λi/n\+λ\)−1\\widehat\{m\}\(\-\\lambda\)=\\tfrac\{1\}\{n\}\\sum\_\{i\}\(\\lambda\_\{i\}/n\+\\lambda\)^\{\-1\}\. The limit statement follows from the weak convergence of the empirical spectral distribution together with the boundedness and continuity ofx↦\(x\+λ\)−1x\\mapsto\(x\+\\lambda\)^\{\-1\}on\[0,∞\)\[0,\\infty\)for fixedλ\>0\\lambda\>0\. ∎

## Appendix BImplementation details

All experiments ran on a single laptop \(one NVIDIA RTX PRO 2000, 8 GB; 16 CPU cores, 64 GB RAM\), in single precision for network training and double precision for all kernel solves and error computations\. The dataset is the distributedStructuralMechanicsfile pair \(40000 samples\); inputs are reduced toℝ41\\mathbb\{R\}^\{41\}after verifying the broadcast structure exactly\. Splits: the training block is samples 1–20000, the test set is samples 20001–40000\. High\-data protocol: a fixed permutation of the training block reserves 1000 samples for validation, 19000 for fitting\. Low\-data protocol: the first 1250 samples of the training block, of which the last 250 form the validation split\. The test set is never touched during development; each configuration is evaluated on it once, after selection on validation\.

#### FNO\.

Width 64, 14 modes per dimension, 4 spectral layers with pointwise linear skips and GELU, lift from 3 channels \(broadcast load, two coordinates\), projection64→128→164\\to 128\\to 1\. AdamW, learning rate1\.5×10−31\.5\\times 10^\{\-3\}\(2×10−32\\times 10^\{\-3\}with batch 256\), weight decay10−610^\{\-6\}, cosine schedule to10−610^\{\-6\}, 300 epochs, batch 256\. 6\.45M parameters\.

#### Transformer \(implemented, not in the reported ensemble\)\.

Encoder: 41 tokens \(load value and 10 sine/cosine pairs of the coordinate\), 5 pre\-norm blocks, dimension 192, 4 heads, MLP ratio 4\. Decoder: grid queries from Fourier features of both coordinates, two cross\-attention blocks, pointwise head192→192→1192\\to 192\\to 1; 3\.16M parameters\. The cross\-attention over 1681 grid queries is the most compute\-intensive of our means, and the hardware available for this study \(a single laptop GPU that proved unstable under sustained load\) did not let us train it to the level of the other members; we therefore exclude it from the reported ensemble and note it as the natural fourth architecture family for a follow\-up on stable hardware\.

#### UNet\.

Three convolutional scales \(41→21→1141\\to 21\\to 11by average pooling, ceil mode\) and a bottleneck at6×66\\times 6; two3×33\\times 3convolutions with GroupNorm\(8\) and SiLU per block; widths48/96/192/38448/96/192/384; bilinear upsampling with skip concatenation;1×11\\times 1output head\. Input channels as for the FNO\. AdamW, learning rate1\.5×10−31\.5\\times 10^\{\-3\}, weight decay10−510^\{\-5\}, cosine schedule, 200 epochs, batch 256\. 4\.38M parameters\.

#### MSE variant\.

The residual MLP trained with mean squared error on standardized targets in place of the metric loss, otherwise identical recipe; it enters the ensemble as a loss\-diversity member and is also the strongest plain MLP we obtain\.

#### MLP\.

Input layer41→102441\\to 1024, three residual SiLU blocks of width 1024, output1024→16811024\\to 1681\. AdamW, learning rate10−310^\{\-3\}, weight decay10−510^\{\-5\}, cosine schedule, 400 epochs, batch 256\. 4\.91M parameters\. A wider variant \(width 1536, five residual blocks, 14\.5M parameters\) is trained under the same recipe\.

#### Refiner\.

The same residual trunk with input layer\(41\+1681\)→1024\(41\+1681\)\\to 1024: the load is concatenated with the kernel method’s stress prediction for that load\. Training uses four\-fold out\-of\-fold kernel predictions for the input channel and full\-data kernel predictions at evaluation; otherwise identical to the MLP\. 6\.64M parameters\. Reflection augmentation flips the load and permutes both the kernel channel and the target by the row\-reversal index\.

#### Common training details\.

Loss: mean over the batch of‖v^−v‖2/‖v‖2\\\|\\hat\{v\}\-v\\\|\_\{2\}/\\\|v\\\|\_\{2\}in original units \(predictions denormalized inside the loss; targets standardized by the per\-pixel training mean and global standard deviation\), except for the MSE variant above\. Reflection augmentation with probability1/21/2; reflection averaging at evaluation\. Model selection: validation error computed every 10 epochs, best checkpoint kept\. A single seed is trained per architecture; the ensemble’s diversity comes from the architectures and losses rather than from reseeding, and the correlation analysis of Section[4\.3](https://arxiv.org/html/2609.00389#S4.SS3)indicates same\-architecture reseeds would be more correlated still\.

#### Stacking\.

Global weights: convex, initialized at the simplex minimizer of the measured second\-moment matrixSS\(Proposition[6\.4](https://arxiv.org/html/2609.00389#S6.Thmtheorem4)\) and polished by a short random search on the validation metric\. Per\-pixel weights: affine in the members at each grid point, ridge parameter10−310^\{\-3\}, fitted on half the validation split and accepted only if they beat the global weights on the held\-out half, then refitted on the full split\.

#### Kernel stages\.

Matérn\-5/25/2on standardized inputs\. Scale grid\{0\.5,1,2,4\}×\\\{0\.5,1,2,4\\\}\\timesthe median pairwise distance \(estimated on 2000 points\), nugget grid\{10−7,10−5,10−3\}\\\{10^\{\-7\},10^\{\-5\},10^\{\-3\}\\\}\(scaled bynn\), tuned on validation using an 8000\-sample subsample of the training set, refit at the chosen pair on the full training set in double precision\. The pure kernel baseline \(Table[2](https://arxiv.org/html/2609.00389#S4.T2)\) uses the grid\{0\.5,0\.75,1,1\.5,2\}×\\\{0\.5,0\.75,1,1\.5,2\\\}\\timesmedian and nuggets\{10−8,10−6,10−4\}\\\{10^\{\-8\},10^\{\-6\},10^\{\-4\}\\\}directly atn=19000n=19000\. Feature stage: penultimate activations of the ensemble members \(FNO: spatially averaged pre\-projection channels, 128; transformer: mean\-pooled encoder tokens, 192; MLP: trunk output, 1024\), concatenated, standardized, same kernel family and tuning\.

#### Uncertainty\.

Posterior standard deviation \([2](https://arxiv.org/html/2609.00389#S6.E2)\) computed from the Cholesky factor ofK\+n​λ​IK\+n\\lambda Iby triangular solves\. Conformal scaling: the 0\.9\-quantile of‖e‖2/Pλ\\\|e\\\|\_\{2\}/P\_\{\\lambda\}on the validation split multipliesPλP\_\{\\lambda\}on test; coverage is the fraction of test samples whose absolute error norm falls below the scaled band\.