Conformal Uncertainty Quantification Guarantees for Neural Operators

arXiv:2608.28515v1 [math.NA] 28 Aug 2026

Conformal Uncertainty Quantification Guarantees for Neural Operators

Tom Stent Department of Mathematics Imperial College London London, SW7 2AZ, UK tom.stent25@imperial.ac.uk

Nicolas Boullé Department of Mathematics Imperial College London London, SW7 2AZ, UK n.boulle@imperial.ac.uk

Abstract Neural operators provide fast surrogate models for approximating operators between function spaces, but their predictions often lack uncertainty quantification. We develop a split conformal framework to guarantee that a calibrated pointwise band around the neural operator output contains the true solution on at least a 1 − γ fraction of the evaluation domain, with probability at least 1 − α over test and calibration inputs, where α, γ ∈ (0, 1). Our method reduces a normalized residual field to its spatial (1 − γ)-quantile and computes a scaling factor using a held-out calibration dataset. We prove marginal coverage guarantees for measurable residual fields defined on arbitrary probability spaces, covering both continuum domains and fixed discretizations. Under mild assumptions on the data distribution, we show that the coverage conditional on the calibration set follows a Beta distribution, which we verify with numerical experiments on Darcy flow and Navier–Stokes equations, where our calibration yields bands consistently tighter than existing corrections while retaining the target coverage.

1

Introduction

Operator learning (Boullé and Townsend, 2024; Kovachki et al., 2023, 2024; Li et al., 2021; Lu et al., 2021) has become an increasingly popular approach within the field of scientific machine learning (Karniadakis et al., 2021) to approximate nonlinear operators, such as solution operators associated with partial differential equations (PDEs), between infinite-dimensional function spaces using pairs of input-output functions. A wide range of neural operator architectures, generalizing neural networks to infinite dimensions, have been recently introduced such as DeepONet (Lu et al., 2021), Fourier neural operator (FNO) (Kovachki et al., 2023; Li et al., 2021), convolutional neural operators (Raonic et al., 2023), U-shaped neural operators (Rahman et al., 2023), multi-resolution architectures (Bouziani and Boullé, 2024; Gupta et al., 2021; He et al., 2024; Li et al., 2020), transformer neural operators (Cao, 2021; Hao et al., 2023), and foundation models (Herde et al., 2024; Wang et al., 2026a). These data-driven PDE surrogate techniques have also been applied to large-scale weather forecasting systems (Bi et al., 2023; Lam et al., 2023; Kurth et al., 2023). While several theoretical works studied approximation and sample complexity properties of neural operators (Boullé et al., 2024, 2023; Kovachki et al., 2021, 2024; Lanthaler et al., 2022), ensuring the reliability of their predictions remains a key challenge. In particular, quantifying uncertainty or estimating prediction errors in operator learning may be crucial for applications like weather forecasting, climate modelling, and engineering design. Several neural operator uncertainty quantification methods have been recently introduced (Psaros et al., 2023; Zou et al., 2024), such as Bayesian neural operators (Magnani et al., 2025a), function-valued Gaussian processes (Magnani et al., 2025b), and probabilistic neural operators trained with proper scoring rules (Bülte et al., 2025). Ensemble methods have also been used to identify high-error regions and improve operator predictions under distribution Preprint.

shift (Mouli et al., 2024). These approaches adapt broader Bayesian and ensemble techniques (Gal and Ghahramani, 2016; Gawlikowski et al., 2023; He et al., 2026; Lakshminarayanan et al., 2017), but the resulting uncertainty estimates usually depend on modelling assumptions, training objectives, or data distribution, and generally lack coverage guarantees. Conformal prediction is a technique that uses nonconformity scores computed on a held-out calibration set to construct a prediction set containing the ground truth solution with high probability (Angelopoulos and Bates, 2023; Angelopoulos et al., 2024; Shafer and Vovk, 2008; Vovk et al., 2005). One computes these scores on calibration inputs and uses an upper empirical quantile as the threshold; this yields finite sample, marginal coverage under exchangeability, i.e., coverage in probability over a new draw from the same data distribution. Recent scientific machine learning methods exploit conformal prediction cellwise, grouping nonconformity scores across pairs of calibration sample and grid location and calibrating a single quantile (Gopakumar et al., 2024, 2025, 2026; Moya et al., 2025; Yu et al., 2026). One limitation of this approach is that the guarantee is marginal over both test solutions and spatial locations, hence controlling the average fraction of covered solution-location pairs, instead of the fraction of the domain covered for each individual solution. Gray et al. (2025); Wang et al. (2026b) construct simultaneous prediction sets for function-valued surrogate outputs to obtain strong spatial guarantees. However, the coverage band must accommodate the worst location on the grid, which may result in uninformative estimates for difficult problems. Then, Millard et al. (2025) extend split conformal prediction to neural operators by calibrating in a discretized function space and lifting the guarantee to the underlying function space, assuming that the discretization map is bilipschitz. Closer to this work is the UQNO framework introduced in (Ma et al., 2024), which trains two neural operators: a predictor for the solution and an estimator for the error field, then calibrates a single scaling factor to guarantee that the ground truth is covered by a prediction band over at least a target fraction of the spatial domain. Ma et al. (2024) employ a Hoeffding correction to lift empirical spatial coverage to the continuum domain. This correction requires i.i.d. spatial evaluation points and therefore does not apply to the deterministic grids used by standard FNO implementations and may overestimate the error. Main contributions. This work rigorously formulates split conformal calibration for neural operators using a generic measurable residual field, modelling the ratio between the prediction error and the estimated one, over an arbitrary probability space. This framework applies to both the continuum and discretized settings, giving finite-sample marginal coverage under exchangeability of the data without concentration corrections. Moreover, under weak assumptions on the data distribution (i.i.d. inputs inducing a continuous distribution of conformal scores), we show that the coverage conditional on the calibration set follows a Beta distribution, which quantifies the variance of the marginal guarantee with respect to the calibration set size. Our numerical experiments on the two-dimensional Darcy flow and Navier–Stokes equations demonstrate the validity of our coverage guarantees.

2

Calibration Set and Coverage Guarantees

Setup. Let D ⊂ Rd be a bounded, Lebesgue measurable, domain, with 0 < Leb(D) < ∞, and define the normalized Lebesgue measure ν(B) := Leb(B)/Leb(D). We consider separable Banach spaces A = A(D; Rda ) and U = U(D; Rdu ), whose elements are measurable vector-valued input and output functions on D. We equip these spaces with their respective Borel sigma algebras B(A) and B(U ), and consider the Borel measurable operators G † , Ĝ : A → U, representing the true PDE solution operator and a pre-trained neural operator (approximating G † ). We study the Borel measurable error estimator Ê : A → U(D; R), which outputs a pointwise scalar estimate of the error of predictions of Ĝ, with Ê(a) ≥ 0 ν-almost everywhere (a.e.) for every a ∈ A. We model Ê as an FNO trained on the quantile loss of the error field of predictions of Ĝ (Ma et al., 2024), however all results apply for any pointwise error estimation technique satisfying the relevant assumptions. Note that FNO architectures consist of compositions of continuous operators, so are themselves continuous (and therefore measurable). Therefore, for Ĝ and Ê represented as such, Borel measurability holds. Moreover, if the underlying PDE problem is Hadamard well-posed between the spaces A and U, then G † : A → U is continuous, hence Borel. We emphasize that continuity is only a sufficient condition as we only require Borel measurability of G † , Ĝ and Ê. 2

Motivation. Let a ∈ A be an input function and consider the neural operator prediction Ĝ(a) ≈ G † (a) ∈ U and the error estimator Ê(a) ≈ ∥G † (a) − Ĝ(a)∥2 . Given x ∈ D, we aim to construct a closed ball prediction band Cλ (a)(x) = B̄(Ĝ(a)(x), λÊ(a)(x)) ⊂ Rdu around the prediction Ĝ(a)(x), with a scaling factor λ ∈ [0, ∞] controlling the radius of the ball. The goal is to choose λ such that the true solution G † (a)(x) lies within this band for a fraction 1 − γ of points x ∈ D (with respect to the measure ν), with probability greater than 1 − α over the draw of a from a distribution µ on A, where α, γ ∈ (0, 1) are fixed. 2.1

Calibration set and scaling factor

Since G † (a), Ĝ(a) and Ê(a) are only defined up to equivalence classes in U, the containment set {x ∈ D : G † (a)(x) ∈ Cλ (a)(x)} is not specified pointwise. We must therefore check both that its ν-measure is unchanged by the choice of representatives, and that this measure is a Borel function of a. It suffices to assume the existence of product-measurable representatives g † , ĝ : A × D → Rdu , ê : A × D → R≥0 of the operators G † , Ĝ and Ê. This assumption holds whenever U embeds continuously in C(D; Rdu ) or in some Lp (D; Rdu ) for 1 ≤ p < ∞, which covers every Lebesgue, Sobolev and Hölder space on D. For Ĝ and Ê this is immediate: FNOs are continuous operators with trigonometric polynomial outputs, so the continuous representative is unique and depends measurably on a. For G † only Borel measurability of the operator is required, which follows from Hadamard well-posedness of the PDE. Definition 2.1 (Prediction band). Let λ ∈ [0, ∞], a ∈ A, and x ∈ D, we define the prediction band Cλ (a)(x) ⊂ Rdu around the prediction Ĝ(a)(x) as  Cλ (a)(x) := y ∈ Rdu : ∥y − Ĝ(a)(x)∥2 ≤ λ Ê(a)(x) = B̄(Ĝ(a)(x), λÊ(a)(x)), with the convention that C∞ (a)(x) := Rdu . The constant λ ∈ [0, ∞] is a “scaling factor”, which uniformly scales the error estimator Ê(a)(x). Lemma A.1 shows that the containment fraction of the domain, i.e., ν({x ∈ D : G † (a)(x) ∈ Cλ (a)(x)}), is independent of the choice of representatives for the operators, so that one can drop the representative notation and write G † , Ĝ, Ê for any fixed choice of product-measurable representatives. Let µ be a probability measure on (A, B(A)), and fix a spatial tolerance γ ∈ (0, 1), the fraction of the domain on which containment may fail for a given input, and a probability tolerance α ∈ (0, 1), the maximum probability that this spatial requirement may fail. Let A1 , . . . , An+1 be exchangeable input functions with common marginal law µ, and denote the calibration set by Dcal = (A1 , . . . , An ). We aim to choose a scaling factor λ̂ = λ̂(Dcal ) such that

PA1 ,...,An+1 ν {x ∈ D : G † (An+1 )(x) ∈ Cλ̂ (An+1 )(x)} ≥ 1 − γ ≥ 1 − α, (1) i.e., the true output for the test input is contained on at least a fraction 1 − γ of D with probability at least 1 − α. This is a marginal guarantee over the joint draw of the calibration set and test input, and exchangeability alone is sufficient. Under the assumption that A1 , . . . , An+1 are independent and identically distributed (i.i.d.), the calibration set is independent of the test input and the marginal probability in Eq. (1) decomposes as h
i EDcal PAn+1 ν({x : G † (An+1 )(x) ∈ Cλ̂ (An+1 )(x)}) ≥ 1 − γ Dcal , (2) an average over calibration sets of the conditional probability of sufficient containment. The marginal guarantee controls this average, not the conditional probability for every realized calibration set. Coverage for arbitrary residual fields. Before specializing to neural operators, Section A.2 establishes the conformal argument in a general setting: given a measurable residual field r : A × X → [0, ∞] on a probability space (X, X , π), we define a scalar score as the smallest threshold that contains at least a 1 − γ fraction of the evaluation space X and calibrate an order statistic of these scores to obtain marginal coverage at least 1 − α. Here, the argument only depends on measurability of the field and nothing else, not the model, the PDE, or even that X be a subset of the domain D. The sections that follow apply Theorem A.8 by specifying the evaluation space, its probability measure, and the relevant residual fields: first for the idealized continuum problem, then for discrete observations. 3

2.2

Continuum coverage guarantee

We first consider the idealized setting in which the truth, prediction, and error estimator can be evaluated at every point x ∈ D, instead of at a finite set of points. To define the scaling factor λ associated with the prediction band Cλ from Definition 2.1, we introduce the continuum residual field ρ : A × D → [0, ∞] as pointwise prediction error normalized by the error estimator: ρ(a, x) =

∥G † (a)(x) − Ĝ(a)(x)∥2 Ê(a)(x)

∈ [0, ∞],

a ∈ A, x ∈ D,

(3)

where ρ(a, x) = 0 if Ê(a)(x) = 0 and G † (a)(x) = Ĝ(a)(x), and ρ(a, x) = ∞ if Ê(a)(x) = 0 and G † (a)(x) ̸= Ĝ(a)(x). Since the operators G † , Ĝ, and Ê are Borel measurable, the residual field ρ is jointly measurable with respect to the product sigma algebra B(A) ⊗ B(D). Then, for every a ∈ A and x ∈ D, G † (a)(x) ∈ Cλ (a)(x) if and only if ρ(a, x) ≤ λ, that is {x ∈ D : ρ(a, x) ≤ λ} = {x ∈ D : G † (a)(x) ∈ Cλ (a)(x)} so that one can apply the results from Section A.2 to obtain a marginal coverage guarantee for the continuum problem. Choice of scaling factor. Given a ∈ A, we respectively define the “containment fraction” wλ (a) of the domain for which the residual is smaller than λ and the corresponding scalar nonconformity score s(a), measuring the spatial (1 − γ)-quantile, as
wλ (a) = ν {x ∈ D : ρ(a, x) ≤ λ} , and s(a) = inf{λ ≥ 0 : wλ (a) ≥ 1 − γ}. (4) Let k = ⌈(1 − α)(n + 1)⌉, and consider calibration inputs A1 , . . . , An ∈ A. The scaling factor λ̂ is chosen as the kth smallest value of the nonconformity scores s(A1 ), . . . , s(An ) if k ≤ n and ∞ otherwise. Lemma A.6 then ensures that the nonconformity score s and the resulting calibration map λ̂ (as a function of the calibration inputs A1 , . . . , An ) are measurable. Remark 2.2 (Ideal error estimator). If Ê(a)(x) = ∥G † (a)(x) − Ĝ(a)(x)∥2 for all x ∈ D, then ρ(a, x) = 1 for every a. Every calibration score then satisfies s(Ai ) = 1, so the calibrated scaling factor is λ̂ = 1, and the prediction band C1 (a)(x) is a ball of radius Ê(a)(x) whose boundary passes exactly through G † (a)(x), achieving coverage over the entire domain D. The following result is the continuum instantiation of the split-conformal guarantee proved in Theorem A.8, with (X, X , π) = (D, B(D), ν) and r = ρ. Theorem 2.3 (Continuum marginal coverage). Suppose that A1 , . . . , An+1 ∼ µ are exchangeable, then

k ≥ 1 − α. PA1 ,...,An+1 ν {x ∈ D : G † (An+1 )(x) ∈ Cλ̂ (An+1 )(x)} ≥ 1 − γ ≥ n+1 Moreover, if the Ai are i.i.d. and the
distribution of s(A1 ) is atomless, then the random variable P wλ̂ (An+1 ) ≥ 1 − γ A1 , . . . , An follows a Beta(k, n + 1 − k) distribution and

k 1 Pµ⊗(n+1) ν {x ∈ D : G † (An+1 )(x) ∈ Cλ̂ (An+1 )(x)} ≥ 1 − γ = ≤1−α+ . n+1 n+1 Remark 2.4 (Finiteness of the scaling factor). If the error estimator is constructed (e.g. with a strictly positive final activation function like softplus Dugas et al. 2000) such that Ê(a)(x) > 0 ν-a.e. for µ-a.e. a ∈ A, then Lemma A.5 implies that the scaling factor satisfies λ̂ < ∞ µ⊗n -a.s. provided ⌈(1 − α)(n + 1)⌉ ≤ n. The guarantee in Theorem 2.3 is a marginal coverage guarantee that is finite-sample and distributionfree. While it requires exchangeability of the calibration sets, it does not require any accuracy assumptions on Ĝ or Ê. The second part of Theorem 2.3 provides a precise characterization of the marginal coverage probability in the i.i.d. case, which is useful for understanding the tightness of the guarantee. This result provides a way to quantify the uncertainty associated with the coverage guarantee and to assess the reliability of the prediction bands generated by the conformal method. The atomless assumption ensures that the nonconformity scores are continuous, which is common in operator learning where the input functions are often drawn from a Gaussian process (Kovachki et al., 2023; Boullé and Townsend, 2024). Finally, we remark that the continuum score in Eq. (4) may not be directly computable from finitely many function evaluations, which motivates the discrete guarantees developed in the next section. 4

2.3

Discrete coverage guarantee

We now consider the practical setting in which the solutions, predictions, and the error estimations are evaluated at a discrete set of points DM = {x1 , . . . , xM } ⊂ D. Rather than measuring containment over the continuum domain, our goal is to guarantee that the solution is contained at a sufficiently large fraction of these grid points. We therefore equip DM with the power-set sigma algebra 2DM PM 1 and the uniform counting probability measure νM (B) = M j=1 δxj (B) for B ⊆ DM . Let ĜM and ÊM denote the prediction and error estimators evaluated on DM . For λ ∈ [0, ∞], we define the  discrete prediction band by CλM (a)(xj ) = y ∈ Rdu : ∥y − ĜM (a)(xj )∥2 ≤ λÊM (a)(xj ) and introduce the discrete residual field ρM : A × DM → [0, ∞) as the pointwise error between the prediction and solution, normalized by the error estimator: ρM (a, xj ) =

∥G † (a)(xj ) − ĜM (a)(xj )∥2 ÊM (a)(xj )

,

a ∈ A, xj ∈ DM .

(5)

Here, we assume that a 7→ G † (a)(xj ) and the coordinate maps of ĜM and ÊM are Borel measurable for every j, and that ÊM (a)(xj ) > 0. Since DM is finite, ρM is jointly measurable with respect to B(A) ⊗ 2DM . Then, for every a ∈ A and xj ∈ DM , G † (a)(xj ) ∈ CλM (a)(xj ) if and only if ρM (a, xj ) ≤ λ, so the results of Section A.2 also apply to the discrete problem. We define the fraction of grid points for which the residual is smaller than λ and the corresponding PM 1 scalar nonconformity score similarly to Section 2.2 as wλM (a) = M j=1 1{ρM (a, xj ) ≤ λ} and M sM (a) = inf{λ ≥ 0 : wλ (a) ≥ 1 − γ}. Unlike the continuum case, this spatial quantile can be computed exactly from finitely many evaluations: if q = ⌈(1−γ)M ⌉, then sM (a) is the qth smallest value among ρM (a, x1 ), . . . , ρM (a, xM ). Let k = ⌈(1 − α)(n + 1)⌉, given calibration inputs A1 , . . . , An , the scaling factor λ̂M is chosen as the kth smallest value of the scores sM (A1 ), . . . , sM (An ) if k ≤ n and ∞ otherwise. The following result is a corollary of the generic split-conformal guarantee proved in Theorem A.8 with (X, X , π) = (DM , 2DM , νM ) and r = ρM . Theorem 2.5 (Discrete marginal coverage). Suppose that A1 , . . . , An+1 ∼ µ are exchangeable, then  
k ≥ 1 − α. P νM {xj ∈ DM : G † (An+1 )(xj ) ∈ Cλ̂M (An+1 )(xj )} ≥ 1 − γ ≥ M n+1 Moreover, if the Ai are i.i.d. and the distribution of sM (A1 ) is atomless, then the random variable P(wλ̂M (An+1 ) ≥ 1 − γ | A1 , . . . , An ) follows a Beta(k, n + 1 − k) distribution and M  
1 k P νM {xj ∈ DM : G † (An+1 )(xj ) ∈ Cλ̂M (An+1 )(xj )} ≥ 1 − γ = ≤1−α+ . M n+1 n+1 Remark 2.6 (Calibrating on the discrete score). While the continuum score s from Section 2.2 may not be computable, Corollary A.9 shows that calibrating on any Borel score σ : A → [0, ∞] satisfying s(a) ≤ σ(a) µ-a.s. still yields marginal coverage at level 1 − α for the continuum residual field ρ. In particular, this opens the door to constructing computable scores that dominate s, thereby providing the continuum guarantee of Theorem 2.3 from a discrete implementation.

3

Function-Valued Split Conformal Calibration

We now translate the coverage guarantees from Section 2 into a concrete calibration procedure. Let Nsample D = {(ai , ui )}i=1 ⊂ A × U be a labelled dataset such that ui = G † (ai ) for 1 ≤ i ≤ Nsample . We randomly split D into three disjoint training, residual, and calibration sets. In practice, we typically choose around 60% of the data for training, 30% for residual, and 10% for calibration, but this depends on the difficulty of the problem. The training and residual datasets are used to train the prediction operator and error estimator Ĝ and Ê, while the calibration set is held out until both operators are frozen. This separation is essential as the finite-sample guarantees from Section 2 treat the operators as fixed and require the calibration and test inputs to be exchangeable. Although splitting reduces the data available to fit each model, operator learning has favorable sample complexity in several settings, including elliptic PDE solution operators and noisy linear operators (Boullé et al., 2023; de Hoop et al., 2023). These results suggest that moderate holdout sets need not substantially degrade predictive accuracy, although the effect remains problem-dependent. 5

Algorithm 1 Function-valued split conformal calibration N

sample Input: Labelled data D = {(ai , ui )}i=1 , spatial tolerance γ ∈ (0, 1), and probability tolerance α ∈ (0, 1) Data splitting and model fitting 1: Split D into three disjoint datasets: Dtrain , Dres , and Dcal = {(ai , ui )}n i=1 2: Train the prediction operator Ĝ to approximate ui (x) = G † (ai )(x) on Dtrain and freeze it 3: Train the error estimator Ê to approximate ∥ui (x) − Ĝ(ai )(x)∥2 on Dres and freeze it Conformal calibration 4: for i = 1, . . . , n do 5: Compute theresidual field ρi (x) = ∥ui (x) − Ĝ(ai )(x)∥2 /Ê(ai )(x) for x ∈ D 6: Set Si = inf λ ≥ 0 : ν({x ∈ D : ρi (x) ≤ λ}) ≥ 1 − γ 7: end for 8: Set k = ⌈(1 − α)(n + 1)⌉ and select λ̂ as the kth smallest value among S1 , . . . , Sn if k ≤ n, and λ̂ = ∞ otherwise 9: return Scaling factor λ̂ and prediction band operator Cλ̂ : (a, x) 7→ B̄(Ĝ(a)(x), λ̂Ê(a)(x))

Once the prediction and error estimators are trained, the calibration procedure in Algorithm 1 computes the residual field ρi for each calibration input ai ∈ Dcal . The scalar score Si is defined as the smallest threshold λ such that the residual field ρi does not exceed λ on at least a fraction 1 − γ of the domain, i.e. the spatial (1 − γ)-quantile of ρi under ν. We then select the scaling factor λ̂ as the kth order statistic of these scores, where k = ⌈(1 − α)(n + 1)⌉. This choice ensures that the resulting prediction band Cλ̂ (An+1 ) satisfies the marginal coverage guarantee in Section 2. While the algorithm is stated for the continuum setting, it applies equally to the discrete setting by replacing D with the discrete set DM and ν with the uniform counting measure νM . In this case, the spatial infimum in line 6 reduces to the qth order statistic of the M residual values, with q = ⌈(1 − γ)M ⌉, and is computed by sorting (see Section 2.3). Given the calibration set, the conditional coverage, i.e., the probability that a test input An+1 is covered, follows a Beta distribution:

Q := PAn+1 ν {x ∈ D : G † (An+1 )(x) ∈ Cλ̂ (An+1 )(x)} ≥ 1 − γ | Dcal ∼ Beta(k, n + 1 − k), under the i.i.d. and atomless-score conditions of Theorems 2.3 and 2.5. This characterization quantifies the variability of the conditional coverage across calibration sets since Var(Q) = k(n+1−k) ∼ α(1−α) . Hence, the standard deviation of the conditional coverage is O(n−1/2 ), (n+1)2 (n+2) n→∞ n which decreases as the calibration-set size increases.

4

Numerical Experiments

We compare our functional conformal calibration (Algorithm 1) with the method of Ma et al. (2024). Both methods use the same retrained operators (ĜM , ÊM ) and calibration data to obtain a fair comparison, and produce pointwise bands ĜM (a)(xj ) ± λ̂ÊM (a)(xj ) for each input a and grid node xj , since in this setting du = 1. The neural operator predictor ĜM and error estimator ÊM are Fourier neural operators (FNOs) of identical architecture, with a final Softplus activation in ÊM to enforce non-negativity of the output. We perform numerical experiments on two standard benchmarks from (Li et al., 2021): a steady 2D Darcy flow (also used in Ma et al. 2024) and a more challenging 2D incompressible Navier–Stokes dataset in vorticity form at ν = 10−4 . Each dataset is split once with a fixed seed into four disjoint subsets as training, residual, calibration, and test sets following Section 3. We train ĜM with relative L2 loss (Kossaifi et al., 2025) and AdamW, then freeze it. We then train ÊM on pointwise residuals ∥G † (a)(xj ) − ĜM (a)(xj )∥2 with pinball loss at level 1 − γ, so that it estimates the pointwise (1 − γ)-quantile of the prediction error. Additional implementation details are available in Section B.1. Given a dataset on a two-dimensional grid at resolution Nmax , we subsample it to obtain coarser grids of resolution N , with M = N 2 number of points. At each resolution, we train a new (ĜM , ÊM ) pair and compare two calibration rules: ours (Algorithm 1) and (Ma et al., 2024), using the same trained operators and calibration data. Throughout, containment fraction denotes the fraction of 6

grid points for which a given input is contained in the prediction band, and coverage is the fraction of test inputs whose containment fraction is greater than 1 − γ. In the experiments that follow, we select γ = α = 0.1; we are interested in a 90% containment of the grid points, over 90% of calibration/test draws. We then perform 3000 random disjoint re-splits of the calibration and test functions, recalibrate each time, and report the empirical distribution of λ̂. The validity of the method is assessed using the same resampling procedure, where we compute the empirical coverage on each re-split. Finally, we evaluate the transfer of a pair calibrated at N < Nmax to Nmax by either spectral interpolation of coarse outputs to fine nodes, or zero-shot super-resolution.

5

Ma et al. Ours

4

100

Density

Calibrated scaling factor ¸^

Two-dimensional Darcy flow. In Fig. 1(left), we report the median calibrated scaling factors along with the 5th–95th percentiles over the 3000 random disjoint re-splits of the calibration and test sets, with the trained operators held fixed. The calibrated factors given by our method are consistently smaller than those of (Ma et al., 2024), across spatial resolutions with ratio ranging from 2.37× at N = 12 to 1.16× at N = 84 (Table p 2). While the correction of (Ma et al., 2024) shrinks with resolution (roughly at a scale of − log α/(2N d )), it remains strictly conservative at practical resolutions. Fig. 1(right) shows that our tighter scaling factors do not compromise validity: the empirical distribution of resampled coverage at N = 60 closely matches the Beta-Binomial law
(BetaBinomial ntest , k, ncal + 1 − k with k = ⌈(1 − α)(ncal + 1)⌉) predicted by Theorem 2.5, indicating near-exact finite-sample calibration on this grid. In contrast, the empirical coverage obtained by Ma et al. (2024) concentrates near 0.98, which ensures validity but is substantially more conservative than the target 1 − α = 0.9.

3 2

80

Beta-binomial (exact) Ours Ma et al.

60 40

1−α

20 1 12 20 30

42

60

0 0.80

84

Resolution N

0.85

0.90

0.95

1.00

Resampled empirical coverage

Figure 1: 2D Darcy flow dataset. Left: median calibrated scaling factor λ̂ with respect to the spatial resolution N over 3000 random disjoint re-splits (ncal = 500, ntest = 1000) along with empirical 5th–95th percentiles; smaller λ̂ indicates tighter bands. Right: empirical distribution of resampled coverage at resolution N = 60 compared with the exact Beta-Binomial law predicted by Theorem 2.5 at the target level 1 − α = 0.9. Our distribution aligns with the predicted law, while (Ma et al., 2024) is shifted right with more conservative and wider bands. The violin plots in Fig. 3 display the same behavior across all spatial resolutions N considered, with complete conservativeness (point mass near 1) at coarse resolutions. Finally, Fig. 4 in Section B.2 evaluates transfer of the coverage guarantees to the fine grid. Under zero-shot super-resolution, we note that the coverage is zero at all calibration resolutions for both methods, while under spectral interpolation, coverage is near zero for N ≤ 42 and increases with the spatial resolution. While the correction of (Ma et al., 2024) is designed to provide exactly this transfer by inflating the scaling factor, it only achieves the desired coverage at resolution N = 84. In this case, the discretization is fine enough so that no correction is needed and our method achieves the target coverage without any inflation. 2D incompressible Navier-Stokes. Given an initial vorticity field, we aim to predict the solution at a later time t = 21. The problem is more challenging than Darcy flow, with a higher prediction error (about 13% relative error) and a more complex spatial structure to capture. We train the same FNO architectures as in Darcy flow, but with a larger number of Fourier modes to capture the finer details of the solution fields. We display in Fig. 2, a visualization of the predicted solution for a test input, along with the pointwise error and the scaled (our calibration method) estimated error field. Additional evaluations in Section B.2 show that the empirical coverage distribution again follows the 7

Input

True solution

True error

1.6

2

Predicted solution

Containment (93.9%)

0

1

Predicted error

0.8

¡2 0

0.0

Figure 2: Navier–Stokes at N = 64: input vorticity, solution at t = 21, FNO prediction, pointwise error, and estimated error field for a calibrated operator pair. The estimated error tracks the spatial error pattern, enabling tight pointwise bands while maintaining the target containment fraction 1 − γ.

Beta-Binomial law predicted by Theorem 2.5 (cf. Figs. 5 and 6), and that our method consistently produces tighter bands than Ma et al. (2024) across all resolutions considered. As observed in Table 3, the ratio of the scaling factors produced by their method (Ma et al., 2024) relative to ours ranges from 1.78× at N = 16 to 1.23× at N = 64, while preserving the target coverage. The empirical coverage distribution for Ma et al. (2024) is again shifted to higher values, confirming its conservative nature.

5

Conclusion

In this paper, we developed a split conformal framework to provide theoretically justified uncertainty quantification guarantees for neural operators. The calibration procedure constructs prediction bands based on a held-out calibration sets that are guaranteed to contain the true solution with a tolerance probability, while allowing for a controlled fraction of spatial failure over the domain. The construction reduces a measurable residual field to a spatial quantile and calibrates a single scaling factor, yielding finite-sample marginal coverage under exchangeability. The proof strategy applies to both continuum domains and fixed grids, with a discrete guarantee that does not rely on random spatial sampling or concentration corrections, and we characterized calibration-set variability via an exact Beta law under i.i.d. sampling and atomless scores. Empirically, on Darcy flow and Navier–Stokes equations, coverage matches the predicted behavior while achieving systematically tighter bands than (Ma et al., 2024), whose corrections are more conservative. The approach is model-agnostic and applies to any neural operator architecture, provided that an appropriate error estimator is available.

Limitations. Our guarantees are marginal over the calibration and test inputs and assume their exchangeability; they do not ensure coverage under distribution shift or guarantee the target coverage conditional on every realized calibration set. Finally, the discrete guarantee is tied to the calibration grid and is not automatically invariant to the spatial resolution. Corollary A.9 offers a framework for lifting it to the continuum, via a computable score dominating the continuum score. Constructing such a score requires regularity of the solution operator and the learned fields, is problem-dependent, and may produce conservative bands. 8

References Anastasios N Angelopoulos and Stephen Bates. Conformal prediction: A gentle introduction. Found. Trends Mach. Learn., 16(4):494–591, 2023. Anastasios N Angelopoulos, Rina Foygel Barber, and Stephen Bates. Theoretical foundations of conformal prediction. arXiv preprint arXiv:2411.11824, 2024. Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Accurate mediumrange global weather forecasting with 3D neural networks. Nature, 619(7970):533–538, 2023. Nicolas Boullé and Alex Townsend. A mathematical guide to operator learning. In Handb. Numer. Anal., volume 25, pages 83–125. Elsevier, 2024. Nicolas Boullé, Diana Halikias, and Alex Townsend. Elliptic PDE learning is provably data-efficient. Proc. Natl. Acad. Sci. USA, 120(39):e2303904120, 2023. Nicolas Boullé, Diana Halikias, Samuel E Otto, and Alex Townsend. Operator learning without the adjoint. J. Mach. Learn. Res., 25(364):1–54, 2024. Nacime Bouziani and Nicolas Boullé. Structure-preserving operator learning. arXiv preprint arXiv:2410.01065, 2024. Christopher Bülte, Philipp Scholl, and Gitta Kutyniok. Probabilistic neural operators for functional uncertainty quantification. Trans. Mach. Learn. Res., 2025. Shuhao Cao. Choose a transformer: Fourier or Galerkin. In Adv. Neural Inf. Process. Syst., volume 34, 2021. Maarten V de Hoop, Nikola B Kovachki, Nicholas H Nelsen, and Andrew M Stuart. Convergence rates for learning linear operators from noisy data. SIAM/ASA J. Uncertain. Quantif., 11(2): 480–513, 2023. Charles Dugas, Yoshua Bengio, François Bélisle, Claude Nadeau, and René Garcia. Incorporating second-order functional knowledge for better option pricing. In Adv. Neural Inf. Process. Syst., volume 13, 2000. Rick Durrett. Probability: Theory and Examples. Cambridge University Press, 5th edition, 2019. Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In Int. Conf. Mach. Learn., pages 1050–1059, 2016. Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, Mohsin Ali, Jongseok Lee, Matthias Humt, Jianxiang Feng, Anna Kruspe, Rudolph Triebel, Peter Jung, Ribana Roscher, Muhammad Shahzad, Wen Yang, Richard Bamler, and Xiao Xiang Zhu. A survey of uncertainty in deep neural networks. Artif. Intell. Rev., 56:1513–1589, 2023. Vignesh Gopakumar, Joel Oskarrson, Ander Gray, Lorenzo Zanisi, Stanislas Pamela, Daniel Giles, Matt Kusner, and Marc Deisenroth. Valid error bars for neural weather models using conformal prediction. arXiv preprint arXiv:2406.14483, 2024. Vignesh Gopakumar, Ander Gray, Lorenzo Zanisi, Timothy Nunn, Daniel Giles, Matt Kusner, Stanislas Pamela, and Marc Peter Deisenroth. Calibrated Physics-Informed Uncertainty Quantification. In Int. Conf. Mach. Learn., 2025. Vignesh Gopakumar, Ander Gray, Joel Oskarsson, Lorenzo Zanisi, Daniel Giles, Matt J Kusner, Stanislas Pamela, and Marc Peter Deisenroth. Uncertainty quantification of surrogate models using conformal prediction. Mach. Learn. Sci. Technol., 7(1):015025, 2026. Ander Gray, Vignesh Gopakumar, Sylvain Rousseau, and Sebastien Destercke. Guaranteed Prediction Sets for Functional Surrogate Models. In Conf. Uncertain. Artif. Intell., pages 1569–1585, 2025. Gaurav Gupta, Xiongye Xiao, and Paul Bogdan. Multiwavelet-based operator learning for differential equations. In Adv. Neural Inf. Process. Syst., volume 34, 2021. 9

Zhongkai Hao, Zhengyi Wang, Hang Su, Chengyang Ying, Yinpeng Dong, Songming Liu, Ze Cheng, Jun Song, and Jun Zhu. GNOT: A general neural operator transformer for operator learning. In Int. Conf. Mach. Learn., 2023. Juncai He, Xinliang Liu, and Jinchao Xu. MgNO: Efficient parameterization of linear operators via multigrid. In Int. Conf. Learn. Represent., volume 2024, pages 53409–53428, 2024. Wenchong He, Zhe Jiang, Tingsong Xiao, Zelin Xu, and Yukun Li. A survey on uncertainty quantification methods for deep learning. ACM Comput. Surv., 58(7):1–35, 2026. Maximilian Herde, Bogdan Raonić, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel De Bezenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for PDEs. In Adv. Neural Inf. Process. Syst., volume 37, pages 72525–72624, 2024. George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nat. Rev. Phys., 3(6):422–440, 2021. Jean Kossaifi, Nikola Kovachki, Zongyi Li, David Pitt, Miguel Liu-Schiaffini, Robert J. George, Boris Bonev, Kamyar Azizzadenesheli, Julius Berner, Valentin Duruisseaux, and Anima Anandkumar. A library for learning neural operators. arXiv preprint arXiv:2412.10354, 2025. Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for Fourier neural operators. J. Mach. Learn. Res., 22(290):1–76, 2021. Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. J. Mach. Learn. Res., 24(89):1–97, 2023. Nikola B Kovachki, Samuel Lanthaler, and Andrew M Stuart. Operator learning: Algorithms and analysis. Handb. Numer. Anal., 25:419–467, 2024. Thorsten Kurth, Shashank Subramanian, Peter Harrington, Jaideep Pathak, Morteza Mardani, David Hall, Andrea Miele, Karthik Kashinath, and Anima Anandkumar. FourCastNet: Accelerating Global High-Resolution Weather Forecasting Using Adaptive Fourier Neural Operators. In Proc. Platf. Adv. Sci. Comput. Conf.. ACM, 2023. Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. In Adv. Neural Inf. Process. Syst., volume 30, 2017. Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, et al. Learning skillful medium-range global weather forecasting. Science, 382(6677):1416–1421, 2023. Samuel Lanthaler, Siddhartha Mishra, and George E Karniadakis. Error estimates for deeponets: A deep learning framework in infinite dimensions. Trans. Math. Appl., 6(1):tnac001, 2022. Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. In Adv. Neural Inf. Process. Syst., volume 33, pages 6755–6766, 2020. Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier Neural Operator for Parametric Partial Differential Equations. In Int. Conf. Learn. Represent., 2021. Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nat. Mach. Intell., 3:218–229, 2021. Ziqi Ma, David Pitt, Kamyar Azizzadenesheli, and Anima Anandkumar. Calibrated Uncertainty Quantification for Operator Learning via Conformal Prediction. Trans. Mach. Learn. Res., 2024. Emilia Magnani, Nicholas Krämer, Runa Eschenhagen, Lorenzo Rosasco, and Philipp Hennig. Approximate Bayesian Neural Operators: Uncertainty Quantification for Parametric PDEs. Trans. Mach. Learn. Res., 2025a. 10

Emilia Magnani, Marvin Pförtner, Tobias Weber, and Philipp Hennig. Linearization turns neural operators into function-valued Gaussian processes. In Int. Conf. Mach. Learn., 2025b. David Millard, Lars Lindemann, and Ali Baheri. Split conformal prediction in the function space with neural operators. arXiv preprint arXiv:2509.04623, 2025. S Chandra Mouli, Danielle C. Maddix, Shima Alizadeh, Gaurav Gupta, Andrew Stuart, Michael W. Mahoney, and Bernie Wang. Using Uncertainty Quantification to Characterize and Improve Out-of-Domain Learning for PDEs. In Int. Conf. Mach. Learn., 2024. Christian Moya, Amirhossein Mollaali, Zecheng Zhang, Lu Lu, and Guang Lin. Conformalizeddeeponet: A distribution-free framework for uncertainty quantification in deep operator networks. Phys. D: Nonlinear Phenom., 471:134418, 2025. Apostolos F Psaros, Xuhui Meng, Zongren Zou, Ling Guo, and George Em Karniadakis. Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons. J. Comput. Phys., 477:111902, 2023. Md Ashiqur Rahman, Zachary E. Ross, and Kamyar Azizzadenesheli. U-NO: U-shaped neural operators. Trans. Mach. Learn. Res., 2023. Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci, Rima Alaifari, Siddhartha Mishra, and Emmanuel De Bézenac. Convolutional neural operators for robust and accurate learning of PDEs. In Adv. Neural Inf. Process. Syst., volume 36, pages 77187–77200, 2023. Glenn Shafer and Vladimir Vovk. A tutorial on conformal prediction. J. Mach. Learn. Res., 9: 371–421, 2008. Vladimir Vovk. Conditional validity of inductive conformal predictors. In Asian Conf. Mach. Learn., pages 475–490. PMLR, 2012. Vladimir Vovk, Alex Gammerman, and Glenn Shafer. Algorithmic Learning in a Random World. Springer, 2005. Hong Wang, Haiyang Xin, Jie Wang, Xuanze Yang, Fei Zha, Yan Jiang, et al. Mixture-of-experts operator transformer for large-scale pde pre-training. In Adv. Neural Inf. Process. Syst., volume 38, pages 31498–31527, 2026a. Weinan Wang, Bowen Gang, and Hao Deng. Operator learning for the 2D incompressible NavierStokes equations: a conformal prediction approach in the data-scarce regime. arXiv preprint arXiv:2606.08654, 2026b. Yifan Yu, Cheuk Hin Ho, and Yangshuai Wang. A conformal prediction framework for uncertainty quantification in physics-informed neural networks. J. Comput. Phys., 561:114979, 2026. Zongren Zou, Xuhui Meng, Apostolos F Psaros, and George E Karniadakis. NeuralUQ: A comprehensive library for uncertainty quantification in neural differential equations and operators. SIAM Rev., 66(1):161–190, 2024.

A

Technical Proofs and Auxiliary Results

A.1

Independence of the representatives

We show that the measure of the containment set defined in Section 2.1 does not depend on the choice of representatives for the operators G † , Ĝ, and Ê. Lemma A.1 (Representative independence). Let (g1† , ĝ1 , ê1 ) and (g2† , ĝ2 , ê2 ) be two families of product-measurable representatives of (G † , Ĝ, Ê), with corresponding prediction bands Cλ1 and Cλ2 . Then for every a ∈ A and λ ∈ [0, ∞],

ν {x ∈ D : g1† (a, x) ∈ Cλ1 (a)(x)} = ν {x ∈ D : g2† (a, x) ∈ Cλ2 (a)(x)} . 11

i Proof. Fix a and λ; the case λ = ∞ is trivial since C∞ (a)(x) = Rdu for i = 1, 2. For λ < ∞ and † i i = 1, 2, set Ei := {x ∈ D : gi (a, x) ∈ Cλ (a)(x)} = {x : ∥gi† (a, x) − ĝi (a, x)∥2 ≤ λ êi (a, x)}. Since each factor is B(A) ⊗ B(D)-measurable, Ei is an a-section of a product-measurable set, hence Ei ∈ B(D) and ν(Ei ) is defined. Because both families represent the same equivalence classes, for each i the equality g1† (a, ·) = g2† (a, ·) holds ν-a.e., and likewise for ĝ and ê. These three a.e. equalities exclude a common ν-null set N ; on D \ N the defining inequalities of E1 and E2 involve identical quantities, so E1 △E2 ⊆ N and ν(E1 ) = ν(E2 ).

A.2

Conformal calibration for a generic residual field

This section provides a generic calibration argument that applies to any residual field over a probability space (X, X , π), including the continuum and discrete residual fields in Section 2. In both cases, there is a nonnegative, jointly measurable field measuring the discrepancy between prediction and truth, a probability measure on the domain over which containment is assessed, a scalar nonconformity score set as a quantile of that field, and a split-conformal calibration procedure that turns exchangeability of the scores into a coverage guarantee. The argument proceeds in four steps. First, we convert a residual field into a scalar score by taking the smallest threshold that contains at least a 1 − γ fraction of the evaluation space. Second, we establish the measurability of this score and show that thresholding it is equivalent to achieving the desired containment fraction. Third, an exchangeablerank argument calibrates the kth order statistic of the scores and yields marginal coverage at least 1 − α. Finally, we show that the same guarantee remains valid when calibration uses a computable score that upper-bounds an ideal score. Let γ, α ∈ (0, 1), and A1 , . . . , An+1 be exchangeable with common law µ on the separable Banach space A, for n ≥ 1. Let r : A × X → [0, ∞] be a nonnegative B(A) ⊗ X -measurable map, which we call a residual field over (X, X , π). The field is evaluated at x ∈ X for each a ∈ A, and measures the discrepancy between a prediction and a ground truth. For a ∈ A and λ ∈ [0, ∞], we say that the field is contained at x if r(a, x) ≤ λ, and denote by
wλr (a) = π {x ∈ X : r(a, x) ≤ λ} (6) the corresponding containment fraction. The quantity wλr (a) measures the fraction (with respect to π) of the evaluation space X for which the residual field is contained by λ. The next lemma ensures that the containment fraction wλr (a) is monotone and right-continuous in λ and measurable in a. Lemma A.2 (Regularity and measurability). For each a ∈ A, the map λ 7→ wλr (a) is nondecreasing r and right-continuous, with w∞ (a) = 1. Moreover, for each λ ∈ [0, ∞], the map a 7→ wλr (a) is Borel measurable. Proof. Let 0 ≤ λ1 ≤ λ2 . Then {x : r(a, x) ≤ λ1 } ⊆ {x : r(a, x) ≤ λ2 }, and monotonicity of π implies that wλr 1 (a) ≤ wλr 2 (a). For the right-continuity, let λ0 ≥ 0 and λm ↓ λ0 . The sets T {x : r(a, x) ≤ λm } are nonincreasing in m, and m∈N {x : r(a, x) ≤ λm } = {x : r(a, x) ≤ λ0 }: the inclusion ⊇ holds since λm ≥ λ0 for all m. Conversely, if x lies in every set then r(a, x) ≤ λm for all m, and letting m → ∞ gives r(a, x) ≤ inf m λm = λ0 . Continuity from above of the probability measure π yields limm→∞ wλr m (a) = wλr 0 (a). Finally, every x ∈ X satisfies r(a, x) ≤ ∞, so r w∞ (a) = π(X) = 1. For the second part, we note that r is B(A) ⊗ X -measurable as a map into [0, ∞], so (a, x) 7→ 1{r(a, x) ≤ λ} is a nonnegative product-measurable function.R As a probability measure, π is σ-finite, so Tonelli’s theorem applies, and the partial integral a 7→ X 1{r(a, x) ≤ λ} dπ(x) = wλr (a) is B(A)-measurable. Given a ∈ A, we define the nonconformity score sr (a) as the infimum of the scaling factors λ for which the containment fraction wλr (a) is greater than 1 − γ: sr (a) = inf{λ ≥ 0 : wλr (a) ≥ 1 − γ},

(7)

where inf ∅ := +∞ by convention. The right-continuity of λ 7→ wλr (a) is crucial and ensures that the infimum in the definition of sr (a) is attained, so that sr (a) ≤ λ if and only if wλr (a) ≥ 1 − γ, as shown by the following lemma. Lemma A.3 (Score-containment equivalence). For every a ∈ A and λ ∈ [0, ∞], sr (a) ≤ λ if and only if wλr (a) ≥ 1 − γ. 12

Proof. If wλr (a) ≥ 1 − γ, then sr (a) ≤ λ by definition (and trivially when λ = ∞). Conversely, suppose sr (a) ≤ λ. If sr (a) = ∞, then λ = ∞ and wλr (a) = 1. Otherwise choose tm ↓ sr (a) with wtrm (a) ≥ 1 − γ. Following Lemma A.2, the right-continuity of λ 7→ wλr (a) implies that wsrr (a) (a) = limm wtrm (a) ≥ 1 − γ; in particular the infimum is attained. Monotonicity then gives wλr (a) ≥ wsrr (a) (a) ≥ 1 − γ. Set k = ⌈(1 − α)(n + 1)⌉ and for an n-tuple v ∈ [0, ∞]n , denote by Tk (v) its kth smallest entry. Split conformal calibration applies this order statistic to the calibration scores. r n Definition A.4 (Calibrated scaling factor). The calibration map ∞] is defined for
λ̂ : A → [0, r r r n (a1 , . . . , an ) ∈ A as λ̂ (a1 , . . . , an ) = Tk s (a1 ), . . . , s (an ) if k ≤ n, and λ̂r (a1 , . . . , an ) = ∞ otherwise.

We abbreviate λ̂r := λ̂r (A1 , . . . , An ) when no confusion arises, and first connect the geometric containment construction to the scalar conformal calibration. Definition A.4 determines the scaling factor provided by the calibration procedure. In the rest of the section, we show that it provides the desired coverage guarantee, i.e., that the residual field is contained on at least a 1 − γ fraction of X with probability at least 1 − α. The convention λ̂r = ∞ guarantees formal coverage when the requested quantile exceeds the calibration sample, but gives a vacuous band. The following mild condition ensures finite scores, and hence a finite calibrated scaling when k ≤ n. Lemma A.5 (Finiteness of the score). Let a ∈ A such that π({x : r(a, x) < ∞}) = 1, then sr (a) < ∞. Moreover, if this holds for all a ∈ A outside a µ-null set and k ≤ n, then λ̂r < ∞ µ⊗n -almost surely. Proof. Fix a ∈ A with π({x : r(a, x) < ∞}) = 1. For m ∈ N, the sets {x : r(a, x) ≤ m}, m ∈ N increase to {x : r(a, x) < ∞}, so by continuity from below of π, we have r wm (a) −−−−→ π({x : r(a, x) < ∞}) = 1 > 1 − γ, m→∞

r where the inequality is strict because γ > 0. Hence, there exists m ∈ N such that wm (a) ≥ 1 − γ, which implies that sr (a) ≤ m < ∞ by Lemma A.3.

For the second statement, suppose the assumption holds for all a ∈ A outside a µ-null set and k ≤ n. Then sr (Ai ) < ∞ a.s. for each i ≤ n by the first part of the lemma. Therefore, maxi≤n sr (Ai ) < ∞
a.s., and since Tk (v) ≤ maxi vi for k ≤ n, we have λ̂r = Tk sr (A1 ), . . . , sr (An ) ≤ maxi≤n sr (Ai ) < ∞ µ⊗n -almost surely. To apply the probabilistic arguments below, we must verify that the score, the calibrated scaling factor, and their comparison event are measurable. We formulate this for a general Borel score σ so that the result also applies to the dominating-score calibration in Corollary A.9. Lemma A.6 (Event measurability). The score sr and the calibration map λ̂r are Borel measurable. More generally, for any Borel score σ : A → [0, ∞], let λ̂σ be the calibration map obtained by replacing sr with σ in Definition A.4. Then λ̂σ is measurable, and  Eσ = (a1 , . . . , an+1 ) ∈ An+1 : sr (an+1 ) ≤ λ̂σ (a1 , . . . , an ) ∈ B(A)⊗(n+1) . Proof. By combining Lemmas A.2 and A.3, we find that the sublevel sets of sr are Borel: {a ∈ A : sr (a) ≤ λ} = {a ∈ A : wλr (a) ≥ 1 − γ} is Borel for every finite λ and equal to A if λ = ∞, hence sr is Borel. Moreover, if k ≤ n, then [ \ {v ∈ [0, ∞]n : Tk (v) ≤ t} = {v ∈ [0, ∞]n : vi ≤ t}. I⊆{1,...,n} i∈I |I|=k

Each set {v ∈ [0, ∞]n : vi ≤ t} is Borel because the coordinate projection v 7→ vi is continuous. Thus the right-hand side is Borel for every t as a finite union of finite intersections of measurable sets, so Tk is measurable. Since σ is Borel, the coordinatewise score map (a1 , . . . , an ) 7→ (σ(a1 ), . . . , σ(an )) is B(A)⊗n -measurable, since a map into a product space is measurable if and only if each coordinate map is, and each coordinate is σ composed with a coordinate 13

projection. Therefore, the composition λ̂σ = Tk ◦ (σ × · · · ×
σ) is measurable. Finally, this implies that the map (a1 , . . . , an+1 ) 7→ sr (an+1 ), λ̂σ (a1 , . . . , an ) ∈ [0, ∞]2 is B(A)⊗(n+1) -measurable, and {(u, v) ∈ [0, ∞]2 : u ≤ v} is closed, hence Borel, in [0, ∞]2 . Then its preimage Eσ lies in B(A)⊗(n+1) . The following lemma formalizes the decomposition of the marginal coverage probability in (1) into an average of conditional coverage probabilities. Conditioning on the calibration set freezes λ̂ at a fixed number, so the only randomness left is the draw of the test input. The conditional probability of coverage is then simply the proportion of inputs whose score falls below that number, that is, the score distribution function evaluated at λ̂. Thus F (λ̂r ) varies across random draws of the calibration set, but is fixed once a particular calibration set is realized. The marginal guarantee controls the mean of this quantity over calibration sets; it does not require the conditional coverage for every realized calibration set to be at least 1 − α. Proposition A.7 (Conditional coverage identity). Suppose (A1 , . . . , An ) is independent of An+1 ∼ µ, and for λ ∈ [0, ∞], denote by F (λ) = µ({a : sr (a) ≤ λ}), the distribution function of the score sr (An+1 ). Then, for almost every realization of the calibration sample (A1 , . . . , An ) with respect to its law (which is µ⊗n when the calibration inputs are i.i.d.),
P wλ̂r r (An+1 ) ≥ 1 − γ A1 , . . . , An = F (λ̂r ),
  hence P sr (An+1 ) ≤ λ̂r = E F (λ̂r ) . Proof. Let V = (A1 , . . . , An ) and define h : An × A → {0, 1} by h(v, w) := 1{sr (w) ≤ λ̂r (v)} for v ∈ An and w ∈ A. This function is bounded and jointly measurable by Lemma A.6. For every (v, w) ∈ An × A, Lemma A.3 applied with λ = λ̂r (v) implies that wλ̂r r (v) (w) ≥ 1 − γ if and only if sr (w) ≤ λ̂r (v). Therefore, we have the following equivalence of events: {wλ̂r r (V ) (An+1 ) ≥ 1 − γ} = {sr (An+1 ) ≤ λ̂r (V )} = {h(V, An+1 ) = 1}. Since V and An+1 are independent and h is bounded and measurable, the standard independence lemma for conditional expectations (Fubini plus a monotone class argument; Durrett 2019, Ex. 4.1.7) implies that the following equality holds almost surely, Z  
E h(V, An+1 ) | V = h(V, w) µ(dw) = µ {w : sr (w) ≤ λ̂r (V )} = F (λ̂r ). A

 The second identity is the tower property of conditional expectation P(Esr ) = E E[1Esr |  A1 , . . . , An ] , where Esr = {sr (An+1 ) ≤ λ̂r } is the event defined in Lemma A.6. We can now apply an exchangeable-rank argument to the scalar scores and translate the result back to the containment score using Lemma A.3. Theorem A.8 (Marginal coverage guarantee). Let r : A × X → [0, ∞] be a residual field over (X, X , π), select γ, α ∈ (0, 1) as the spatial and probability tolerances, and suppose A1 , . . . , An+1 are exchangeable. Let k = ⌈(n + 1)(1 − α)⌉ and consider the scaling factor λ̂r defined as the kth smallest value among the scores sr (A1 ), . . . , sr (An ) if k ≤ n, and λ̂r = ∞ otherwise. Then,
P wλ̂r r (An+1 ) ≥ 1 − γ ≥

k ≥ 1 − α, n+1

r that is the coverage is greater than 1 − α. Moreover, if the Ai are i.i.d. and
the law of s (A1 ) is r atomless, then the random variable P wλ̂r (An+1 ) ≥ 1 − γ A1 , . . . , An follows a Beta(k, n + 1 − k) distribution, in particular

P wλ̂r r (An+1 ) ≥ 1 − γ = E[P wλ̂r r (An+1 ) ≥ 1 − γ 14

A1 , . . . , A n ] =

k 1 ≤ 1−α+ . n+1 n+1

Proof. Following Lemma A.6 with σ = sr , the set Esr defined in the lemma’s statement is an event in B(A)⊗(n+1) which coincides with {wλ̂r r (An+1 ) ≥ 1 − γ} by Lemma A.3. We define Si := sr (Ai ) so that Esr = {Sn+1 ≤ λ̂r }. If k > n then λ̂r = ∞ and P(Esr ) = 1 trivially. Otherwise, λ̂r = Tk (S1 , . . . , Sn ), the kth smallest value of the calibration scores. The variables S1 , . . . , Sn+1 are exchangeable, being the image of an exchangeable vector under the fixed Borel map sr applied coordinatewise. The split conformal coverage theorem (Angelopoulos et al., 2024, Fact 2.15(i) & Prop. 3.12) implies that k ≥ 1 − α. n+1 For the second claim, if the Ai are i.i.d. then S1 , . . . , Sn+1 are i.i.d., and since their law is atomless they are pairwise distinct almost surely and the probability distribution function, defined in Propo
sition A.7, F (λ̂r ) = P wλ̂r r (An+1 ) ≥ 1 − γ A1 , . . . , An is continuous. Then, (Angelopoulos

P {wλ̂r r (An+1 ) ≥ 1 − γ} = P Sn+1 ≤ Tk (S1 , . . . , Sn ) ≥

et al., 2024, Sec. 4.1) implies that F (λ̂r ) ∼ Beta(k, n + 1 − k). Theorem A.8 quantifies the calibration-set to calibration-set variability that the marginal guarantee averages over, and is the starting point for calibration-conditional guarantees (Vovk, 2012). In certain applications, the idealized residual field may induce a nonconformity score that cannot be computed directly. However, if this score admits a computable upper bound, the following corollary shows that calibrating on the upper bound preserves the desired coverage guarantee, although the resulting prediction band may be conservative. Corollary A.9 (Calibration with a dominating score). Let r : A × X → [0, ∞] be a residual field over (X, X , π) and σ : A → [0, ∞] be Borel with sr (a) ≤ σ(a) for µ-a.e. a ∈ A. Calibrate on σ, i.e., set λ̂σ as the kth smallest value among σ(A1 ), . . . , σ(An ) if k ≤ n and λ̂σ := ∞ otherwise. Then if A1 , . . . , An+1 are exchangeable,
P wλ̂r σ (An+1 ) ≥ 1 − γ ≥ 1 − α. Proof. The variables σ(A1 ), . . . , σ(An+1 ) are exchangeable, being the coordinatewise image of an exchangeable vector under the fixed Borel map σ, so exactly as in Theorem A.8, P σ(An+1 ) ≤
λ̂σ ≥ k/(n + 1) ≥ 1 − α. Up to a µ⊗(n+1) -null set, on this event sr (An+1 ) ≤ σ(An+1 ) ≤ λ̂σ , and Lemma A.3 implies that wλ̂r σ (An+1 ) ≥ 1 − γ. Measurability of the events follows from Lemma A.6. Remark A.10 (One-sidedness of the calibration). If a computable quantity is guaranteed to overestimate the true score of µ-almost everywhere input, then calibrating on it inherits the coverage guarantee for the true residual field. No upper bound accompanies Corollary A.9, since the domination sr ≤ σ is one-sided, so the coverage may exceed 1 − α by more than 1/(n + 1). This conservativeness is the price of calibrating on an alternative score. The results of this section depend only on the measurable structure of A, the fact that π is a probability measure on (X, X ), and the joint measurability of the residual field. They are otherwise agnostic to the models and error estimators used to construct that field, and therefore apply equally to non-FNO architectures, alternative pointwise error estimators, non-uniform grids, and evaluation spaces unrelated to the physical domain D. These results guarantee coverage but do not control the magnitude of λ̂r or the width of the resulting uncertainty bands. Such efficiency depends on how informative the chosen residual field is, rather than on the conformal calibration argument itself.

B

Implementation and Additional Results

B.1

Implementation Details

Data. Both problems use the benchmark datasets released with the FNO (Li et al., 2021). For Darcy flow we use the 421×421 steady-state data (the same data used by Ma et al., 2024), with constant forcing f ≡ 1 and zero Dirichlet boundary conditions. We discard the duplicated final row and column to obtain Nmax = 420 points per axis. For Navier-Stokes we use the ν = 10−4 vorticity data at 64 × 64 (10 000 trajectories), together with the companion set of 20 trajectories at 256 × 256 used as the fine 15

reference grid, and learn the map from vorticity at t = 11 to t = 21. Working resolutions are obtained by strided subsampling of the reference grid, so every resolution sees the same functions. Each dataset is partitioned once, with a fixed seed (42), into four disjoint subsets (Table 1). The partition is drawn uniformly at random over sample indices. The Navier–Stokes data are from the FNO release (https://drive.google.com/drive/folders/1UnbQh2WWc6knEHbLn-ZaXrKUZhp7pjt-, accessed August 2026), the Darcy data were obtained via the neuraloperator library’s DarcyDataset. Architecture and training. The prediction and error operators share the architecture of Table 1 and differ only in their output activation: the error operator applies a final Softplus to ensure nonnegativity. The prediction operator is trained with the relative L2 loss from the Neural Operator library (Kossaifi et al., 2025), then frozen. The error operator is trained on its pointwise residual fields with the pinball loss at level 1 − γ = 0.9. Both use AdamW with a constant learning rate and no schedule. Error targets are rescaled by a single constant computed from the error training set before fitting, for numerical conditioning. The constant is fixed before calibration and absorbed into the calibrated factor, so it does not affect validity. Hyperparameters were not tuned per method or per calibration rule: a single configuration was chosen per problem so that all comparisons are between calibration rules applied to identical trained operators. Mode counts and epochs were chosen following standard FNO settings. Table 1: Configuration for both problems. Fourier modes are given per axis and are capped by the grid size at each resolution. Darcy flow Navier–Stokes Resolutions N Fourier modes Hidden channels Fourier layers Lifting / projection Optimiser Learning rate Weight decay Batch size Epochs (γ, α)

12, 20, 30, 42, 60, 84 10 at N =12; 16 otherwise 32 4 64 / 64 AdamW 10−3 10−4 16 75 (0.1, 0.1)

16, 32, 64 16, 24, 32 32 4 64 / 64 AdamW 10−3 10−4 32 75 (0.1, 0.1)

Prediction training Error training Calibration Test Resampling split

2000 1000 500 1000 (500, 1000)

5000 3000 1000 1000 (500, 1500)

Calibration and evaluation. Both calibration rules are applied to the same trained operators and the same calibration set, so all reported differences are attributable to the rules alone. For the resampled validity and efficiency experiments, the calibration and test sets are pooled and re-split T = 3000 times at the sizes given in Table 1, using identical permutations for both methods. Compute. All training and evaluation was run on CPU on a single 2024 MacBook Pro (Apple M4 chip, 24GB memory). No GPU or cluster resources were used. Training one prediction/error pair takes approximately 5 minutes at the coarsest and 90 minutes at the finest resolution. Calibration and the T = 3000 resampling experiment take under 3 minutes in total. The fine resolution transfer test of Darcy flow, however took approximately 90 minutes to complete. B.2

Additional Numerical Results

This appendix complements Section 4 with per-resolution tables and distributional diagnostics. As in the main text, N denotes resolution per spatial direction and M = N 2 denotes the number of grid points in 2D. On each fixed grid, the relevant validity statement is Theorem 2.5. 16

Table 2: Darcy flow: calibrated factor λ̂ (single calibration split), mean bandwidth B, and width ratio Pntest 1 PM 1 between methods at each resolution. Mean bandwidth is B := ntest i=1 M j=1 λ̂ ÊM (Ai )(xj ). Ours N 12 20 30 42 60 84

λ̂ 1.84 1.59 1.48 1.64 1.42 1.54

B 9.43 × 10−4 6.90 × 10−4 5.52 × 10−4 4.62 × 10−4 5.41 × 10−4 5.04 × 10−4

λ̂ 4.38 2.24 1.92 2.05 1.66 1.79

Ma et al. B 2.24 × 10−3 9.71 × 10−4 7.13 × 10−4 5.77 × 10−4 6.34 × 10−4 5.87 × 10−4

ratio 2.37 1.41 1.29 1.25 1.17 1.16

Resampled empirical coverage

Table 2 quantifies the efficiency gap from the main text on a single split. Our bands are uniformly tighter, from 2.37× at N = 12 to about 1.16–1.17× at the finest resolutions. The correction term used by Ma et al. (2024) decreases with resolution but remains conservative at practical N . The mean bandwidth B is not monotone in N because operators are retrained independently at each resolution, so model quality and discretization both affect this quantity. By contrast, the per-resolution width ratio is smoother because both methods are applied to the same trained operators and residual fields. 1.000 0.975 0.950 0.925 0.900 0.875 0.850

Ours Ma et al. 1¡®

0.825 12

20

30

42

60

84

Resolution N Figure 3: Darcy flow: distribution of re-sampled coverage by resolution and method, with means and target 1 − α. Our distributions are centred near the target at each N , consistent with Theorem 2.5; those of Ma et al. (2024) are shifted upward and collapse to coverage 1 at the coarsest resolution.

Table 3: Navier-Stokes: calibrated factor λ̂, mean bandwidth B, and the width ratio between the methods, per resolution. Ours Ma et al. N λ̂ B λ̂ B ratio 16 1.92 0.419 3.42 0.747 1.78 32 1.76 0.354 2.34 0.472 1.33 64 1.72 0.354 2.11 0.436 1.23 Across both violin plots (Figs. 3 and 6), our method is calibrated around the target across resolutions, while (Ma et al., 2024) is consistently above target. This is the expected efficiency-validity trade-off: both methods are valid, but the comparator is conservative. 17

Zero-shot super-resolution High-res empirical coverage

High-res empirical coverage

Spectral interpolation 1.0 0.8 0.6

Ma et al. Ours

0.4 0.2 0.0 12 20

30

42

60

1.0 0.8 0.6

Ma et al. Ours

0.4 0.2 0.0

84

12 20

Resolution N

30

42

60

84

Resolution N

4

Ma et al. Ours

200

Density

Calibrated scaling factor ¸^

Figure 4: Darcy flow: coverage on the fine reference grid (Nmax = 420) versus calibration resolution. Left: exact spectral interpolation of calibrated coarse-grid bands. Right: zero-shot super-resolution. This evaluates cross-grid transfer, which is outside the scope of Theorem 2.5 (a fixed-grid guarantee). The correction in Ma et al. (2024) increases width (up to 2.4×) yet reaches the target only at N = 84, and zero-shot coverage is 0 at all resolutions.

3

2

150

Beta-binomial (exact) Ours Ma et al.

100 50

1 16

32

0 0.80

64

Resolution N

0.85

0.90

0.95

1.00

Resampled empirical coverage

Figure 5: Navier–Stokes. Left: resampled median calibrated scaling factor λ̂ versus resolution N , with empirical 5th–95th percentile bars. Right: resampled coverage at N = 64 against the Beta-Binomial reference implied by Theorem 2.5 and the target 1 − α. For Navier–Stokes, Fig. 5 and Table 3 mirror the Darcy conclusions. Efficiency improves at every resolution (ratios 1.78×, 1.33×, 1.23×), and the re-sampled coverage for our method tracks the theorem-predicted Beta-Binomial reference closely. The distribution for (Ma et al., 2024) is shifted to higher coverage, confirming conservative validity with wider bands.

18

Resampled empirical coverage

1.00 0.98 0.96 0.94 0.92 0.90 0.88 Ours Ma et al. 1¡®

0.86 0.84 16

32

64

Resolution N Figure 6: Navier–Stokes: distribution of re-sampled coverage by resolution and method, with means and target 1 − α. Our distributions centre near the target, while those of Ma et al. (2024) are shifted upward with smaller spread, indicating conservative validity.

19

Related documents

• Conformal Uncertainty Quantification Guarantees for Neural Operators#792516
• Conformal Uncertainty Quantification Guarantees for Neural Operators#788516
• Hierarchical Group-Conditional Conformal Risk Control for Selective Prediction in Language Models#405700
• Topological Neural Operators#267636
• Deep Q-Learning on Hölder Spaces#280216