That's a brilliant suggestion. We have crafted a more intuitive explanation here.

Our framework makes uncertainty quantification for neural-PDE solvers more intuitive and practically useful by leveraging two key insights. First, we use the PDE residual (how well the solution satisfies the underlying physics equations) as a measure of uncertainty, rather than traditional error metrics that require ground truth data. This means we can evaluate the quality of predictions without needing expensive simulation data - if a predicted solution violates conservation laws by a large amount, we know it's likely to be inaccurate. Second, we calibrate these physics-based error estimates using conformal prediction, which provides statistical guarantees about the uncertainty bounds. This is like having a physics-aware "confidence score" that tells us how much we can trust each prediction.

The practical value becomes clear in applications - imagine using a neural-PDE solver for weather forecasting or fusion reactor control. Our method provides calibrated error bounds that tell you when predictions are physically implausible and should be double-checked with more expensive traditional methods. We offer two complementary ways to quantify uncertainty: "marginal" coverage that gives error bars for each point in space and time (useful for identifying specific problematic regions), and "joint" coverage that evaluates entire predictions (useful for filtering out globally unreliable solutions). This makes the uncertainty estimates both theoretically rigorous and practically actionable for scientists and engineers using these models.

We thank the reviewer for this interesting question and for giving us the opportunity to elaborate further on an interesting point.

Though our focus within this paper has been on quantifying the misalignment of the model with the equations in the domain of the problem (as mentioned in section 4.3), we should account for the initial and boundary conditions.

A well-defined PDE is characterised by the equations on the domain, the initial condition across the domain at $t=0$ and the boundary conditions, reflecting the physics at the boundary. Within a neural-PDE setting, the initial condition does not need to be enforced or measured for as the neural-PDE is set up as an initial-value problem, taking in the initial state to autoregressively evolve the later time-steps and hence does not come under the purview of the neural-PDE's outputs. The boundary conditions, whether Dirichlet, Neumann or periodic, follow a residual structure as outlined in eqn 2 within the paper, allowing us to use it as a PRE-like nonconformity score for performing conformal prediction. In all the problems we have under consideration, the PDEs are modelled under periodic boundary conditions:

$\frac{\partial u }{\partial X} = 0 ; X \in \partial \Omega$

By deploying the above equation as the PRE across the boundary, we can obtain error bars over the boundary conditions as well. Within figure 21 of the revised paper, we demonstrate the error bars obtained by using the boundary conditions as the PRE nonconformity scores for the Navier-Stokes equations.

We have included the treatment of boundary conditions as a section within the appendix (Appendix J).

We thank the reviewer for the helpful comments and the opportunity for further discussion. As we are approaching the last day allowed for PDF revision, could you kindly let us know if there are any other concerns that we should address?

We would like to thank the reviewer for this question as it puts emphasis on an important aspect of neural-PDEs and our method.

As demonstrated in [1], the discretisation of the inputs and hence model outputs plays an important role in the accuracy of the neural-PDE solvers. Though the neural operators are constructed for discretisation-invariant behaviour due to the band-limited nature of the functions, they often exhibit discretisation-convergent behaviour rather than be fully discretisation-invariant. This is of particular importance in the temporal dimensions as these neural-PDE models utilise a discrete, autoregressive-based time-stepping and is baked into the model within its training regime [2]. Due to lack of control in the discretisation within the temporal domain $(dt)$, the PRE estimates tend to have higher numerical errors as well. In figure 9 of the revised paper, we visualise the evaluation of finite difference in 2D+time as a 3D convolution. The finite difference stencil i.e. the convolutional kernel has a unit discretisation of $dx$, $dy$ and $dt$ associated with the problem and is applied over the signal i.e. the output from the neural-PDE $u$ spanning the domain $x,y,t$, where $x \in [0, X], \\; y \in [0, Y], \\; t \in [0, T]$.

We thank the reviewer for raising this discussion and believe that this represents an important component of using our method and have added a subsection explaining the role of discretisation in the appendix D.

[1] F. Bartolucci, E. de Bezenac, B. Raonic, R. Molinaro, S. Mishra, and R. Alaifari, ``Representation Equivalent Neural Operators: a Framework for Alias-free Operator Learning,'' in Proceedings of the 37th Conference on Neural Information Processing Systems (NeurIPS), 2023

[2] M. McCabe, P. Harrington, S. Subramanian, and J. Brown, ``Towards Stability of Autoregressive Neural Operators,'' Transactions on Machine Learning Research, 2023

We would like to humbly disagree with this statement as we believe that there is a novelty in the nonconformity score and how that extends into performing data-free CP. As reviewer mcYW points out, the paper touches on important topics within the SciML community and we believe that this work leads to practical applications of UQ and active learning. The PRE allows to obtaining calibrated uncertainties for neural operator models, while current UQ methods for such models do not provide guaranteed coverage. The calibration is data-free, which addresses the well-known usual limitation of conformal prediction, which for other nonconformity scores relies on potentially expensive calibration data. In the table comparing various features of different UQ methods (Table 1 in Appendix C of the revised paper), the novelty of our work is further highlighted.

Wave Equation - Coverage measured for 2σ(∼ 95%)

Navier-Stokes Equations - Coverage measured for 2σ(∼ 95%)

MHD Equations - Coverage measured for 2σ(∼ 95%)

Within this experiment, we explore the modelling of plasma for Nuclear Fusion, where surrogate models are being used to determine the real-time stability of the reactor. These surrogate models mimic complex numerical codes (such as JOREK / STORM) and predict whether the reactor would enter into a catastrophic disruption event. Hence, being able to provide inexpensive and immediate yet valid error bars is of utmost urgency as it allows you to isolate events and prevent significant damage to the reactor.

In [1], the authors model the evolution of plasma blobs within a Fusion reactor (known as a Tokamak) following reduced magnetohydrodynamics. Plasma, characterised by density $\rho$, electric potential $\phi$ and Temperature $T$ under the absence of magnetic pressure confining it in place, moves radially outward to the wall of the reactor driven by its kinetic pressure. We demonstrate the ability to scale our method by applying it to obtain valid error bars across the multi-variable FNO trained for plasma modelling.

As shown in figure 23 of the revised paper, our method can capture the model error across a range of predictions and can devise error bars that provide guaranteed coverage without needing any additional data. In figure 23(a), we demonstrate the absolute error in the model prediction of the temperature evolution, and correlate that with the PRE over the temperature equations in figure 23(b). The details of the new experiment has been added to Appendix L of the revised paper.

[1] V. Gopakumar, S. Pamela, L. Zanisi, Z. Li, A. Gray, D. Brennand, N. Bhatia, G. Stathopoulos, M. Kusner, M. P. Deisenroth, A. Anandkumar, the JOREK Team, and MAST Team, "Plasma surrogate modelling using Fourier neural operators," Nuclear Fusion, vol. 64, no. 5, p. 056025, Apr. 2024, doi: 10.1088/1741-4326/ad313a

We did not include comparisons to other UQ methods as our method involved quantifying the physics residual error and obtaining its bounds as opposed to directly quantifying the model's predictive error. However, considering the concerns raised by the reviewer regarding providing baseline comparisons and the effectiveness such a comparison provides to our claims, we have conducted experiments to demonstrate how our method stands out as being data-free, not requiring any modifications or sampling, and helps obtain guaranteed coverage bounds in a physics-informed manner. We have added them to the revised paper under Appendix C.

We thank the reviewer for the helpful comments and the opportunity for further discussion. As we are approaching the last day allowed for PDF revision, could you kindly let us know if there are any other concerns that we should address?

Q: Both Figure 4 and Figure 5 illustrate very large confidence intervals for joint CP which would not correspond to 95% coverage. For marginal CP, it is hard to visually tell what coverage is. The same holds for the later figures. The paper should provide additional metrics to evaluate the results quantitatively.

A: We thank the reviewer for this comment. Visualising coverage for marginal CP is difficult as it is done individually for each cell in the spatio-temporal output tensor. In figures 4 and 5 it is shown for a single example at a fixed point in time. Figure 2 checks whether coverage is valid for each error cell and averages across the domain. To provide further clarity, while comparing against other UQ measures, we showcase the coverage obtained for $\alpha=0.05$, averaged across different samples of the inputs along with the l2 fit of the model, both within and out of distribution cases while comparing against other UQ methods in Tables 2-4 within the appendix section C.

Though we disagree with the reviewer's comments that our work only touches on toy problems, we are thrilled at the opportunity to demonstrate further experiments of our method to a more complex scenario.

In [1], the authors model the evolution of plasma blobs within a Fusion reactor (known as a Tokamak) following reduced magnetohydrodynamics. Plasma, characterised by density $\rho$, electric potential $\phi$ and Temperature $T$ under the absence of magnetic pressure confining it in place, moves radially outward to the wall of the reactor driven by its kinetic pressure. We demonstrate the ability to scale our method by applying it to obtain valid error bars across the multi-variable FNO trained for plasma modelling.

As shown in figure 23 of the revised paper, our method can capture the model error across a range of predictions and can devise error bars that provide guaranteed coverage without needing any additional data. In figure 23(a), we demonstrate the absolute error in the model prediction of the temperature evolution, and correlate that with the PRE over the temperature equations in figure 23(b). The details of the new experiment has been added to Appendix L of the revised paper.

[1] V. Gopakumar, S. Pamela, L. Zanisi, Z. Li, A. Gray, D. Brennand, N. Bhatia, G. Stathopoulos, M. Kusner, M. P. Deisenroth, A. Anandkumar, the JOREK Team, and MAST Team, "Plasma surrogate modelling using Fourier neural operators," Nuclear Fusion, vol. 64, no. 5, p. 056025, Apr. 2024, doi: 10.1088/1741-4326/ad313a

We did not include comparisons to other UQ methods as our method involved quantifying the physics residual error and obtaining its bounds as opposed to directly quantifying the model's predictive error. However, considering the concerns raised by the reviewer regarding providing baseline comparisons and the effectiveness such a comparison provides to our claims, we have conducted experiments to demonstrate how our method stands out as being data-free, not requiring any modifications or sampling, and helps obtain guaranteed coverage bounds in a physics-informed manner. We have added them to the revised paper under Appendix C.

Q: Figure 4 and 5 provide a single illustrative example and not a proper Monte Carlo study to illustrate the calibration and coverage of the method. Considering that the paper suggests the method is computationally cheap, there is no reason to not provide more extensive results supporting the paper's claims.

A: Thank you for this comment, but we believe there is a slight misunderstanding. Although the reviewer is correct that Figures 4. and 5. show single examples, a detailed Monte Carlo validation of our method has been performed. It is what Figure 3 of the main text shows. This figure shows how the empirical coverage, the true coverage computed using Monte Carlo simulation, changes as we change the $\alpha$-level of our constructed confidence intervals. This is an empirical / Monte Carlo validation of equation (5) of the main paper (marginal coverage of conformal prediction). While comparing against other UQ methods in tables 2-4, we conduct further coverage checks by performing MC sampling across a larger dataset.

Q: On line 155 there is a claim: “We introduce a novel data-free nonconformity score for Conformal Prediction (CP) to obtain statistically valid and guaranteed error bounds for neural PDE solvers.” Theoretical support should be provided to this claim.

A: We thank the reviewer for this comment and at this opportunity we would like to highlight that our work extends from the principles of CP, where we extend the coverage guarantees to the residual space. Our work reformulates the prediction sets to be input-independent, allowing us to perform data-free CP. We had explored this by mathematical construction in appendix A, but valuing the feedback, we have formulated this into a Theorem. We have added this to the appendix of our paper.

Preliminaries:

Let $D: \mathbb{R}^m \rightarrow \mathbb{R}^m$ be a physics residual operator mapping a function to its PDE residual value, where: $\\{X_i\\}_{i=1}^n$ is the calibration set, $\hat{f}$ is the model, $\hat{q}^\alpha$ is estimated as the $\lceil(n+1)(1-\alpha)\rceil/n$ -quantile of

$\\{|D(\hat{f}(X_i))|\\}_{i=1}^n$

Theorem:

If the residuals $\\{D(\hat{f}(X_i))\\}_{i=1}^{n+1}$

are exchangeable random variables, then for any significance level $\alpha \in (0,1)$ and any new input $X{n+1}$ we have the following coverage guarantee: $\mathbb{P}(|\mathcal{D}(\hat{f}(X_{n+1}))| \in C_\alpha) \geq 1 - \alpha , \\; C_\alpha = [-\hat{q}\alpha, \hat{q}\alpha]$

Proof:

Let $R_i = |D(\hat{f}(X_i))|$ for $i=1,\ldots,n+1$. We have, by assumption, $(R_1,\ldots,R_n,R_{n+1})$ is an exchangeable sequence. Define the rank $\pi$ of $R_{n+1}$ w.r.t. all other residuals: $\pi(R_{n+1}) = |{i=1,\ldots,n+1: R_i \leq R_{n+1}}|$ By exchangeability, the rank $\pi(R_{n+1})$ is uniformly distributed over ${1,\ldots,n+1}$. Therefore, $P(\pi(R_{n+1}) \leq \lceil(n+1)(1-\alpha)\rceil) = \frac{\lceil(n+1)(1-\alpha)\rceil}{n} \geq 1-\alpha.$ By construction of $\hat{q}^\alpha$ we have that, ${\pi(R_{n+1}) \leq \lceil(n+1)(1-\alpha)\rceil} \subseteq {R_{n+1} \leq \hat{q}^\alpha}.$ Putting this together, $P(|D(\hat{f}(X_{n+1}))| \leq \hat{q}^\alpha) = P(R_{n+1} \leq \hat{q}^\alpha) \geq 1-\alpha,$ which completes the proof. ■

Q: The paper does not distinguish between epistemic and aleatoric uncertainties and should be more explicit and precise about what type of uncertainty quantification it provides.

A: We thank the reviewer for this question, it is the kind of uncertainty that conformal prediction characterises. From one perspective, conformal prediction can be viewed as aleatoric uncertainty, since we construct confidence intervals with respect to a specific probability distribution (the distribution the calibration data is drawn from). However, one could similarly argue the uncertainty is epistemic, since we are model the neural network's error, and use confidence intervals to-do so (used for scalar but unknown quantities, and therefore can be seen as Epistemic). The author's feeling is the second.

The characterisation may not be useful with regards to this paper. As far as we understand, the characterisation is most useful when both types of uncertainty is present in the same problem (so that they don't become confounded, and that appropriate methods are used [1]). Perhaps this could be elaborated in the main text, however since this paper is not on fundamental statistics or fundamental Uncertainty Quantification, we feel the discussion is perhaps slightly out of scope.

[1] S. Ferson, and L. Ginzburg, "Different methods are needed to propagate ignorance and variability", Reliability Engineering & System Safety, Elsevier, 1996, https://doi.org/10.1016/S0951-8320(96)00071-3

Q: “The majority of literature in UQ for neural PDE solvers has been looking at Bayesian methods, such as dropout, Bayesian neural networks and Monte Carlo methods (Geneva & Zabaras, 2020; Zou et al., 2024; Psaros et al., 2023), which lack guarantees or are computationally expensive.”

A: Thank you for the opportunity to clarify our statement regarding uncertainty quantification methods in neural PDE solvers. The primary methods - MC Dropout, Bayesian Neural Networks, and Deep Ensembles - each present significant practical challenges despite their theoretical appeal. MC Dropout requires numerous forward passes for reliable uncertainty estimates, which introduces substantial computational overhead during inference. Bayesian Neural Networks, while theoretically well-founded, demand extensive computational resources for MCMC or variational inference sampling and significantly increase memory requirements due to maintaining weight distributions. Deep Ensembles, though conceptually straightforward, necessitate training and storing multiple complete models, with costs scaling linearly with ensemble size. All these methods, being approximate Bayesian methods share common limitations: they lack theoretical guarantees on uncertainty coverage (especially for out-of-distribution samples), require extensive sampling and modifications to the model or training regime. These challenges become particularly acute in real-world applications where both computational efficiency and reliable uncertainty estimates are crucial. We appreciate the reviewer's attention to this point and would be happy to provide additional technical details if needed.

Wave Equation - Coverage measured for 2σ(∼ 95%)

Navier-Stokes Equations - Coverage measured for 2σ(∼ 95%)

MHD Equations - Coverage measured for 2σ(∼ 95%)

We thank the reviewer for the helpful comments and the opportunity for further discussion. As we are approaching the last day allowed for PDF revision, could you kindly let us know if there are any other concerns that we should address?

We thank the reviewer for this question, which touches upon an important aspect of neural PDEs: their fit not just to the data but to the underlying physics. The reviewer is spot-on in assuming that PRE-CP will still provide guaranteed coverage with considerably wider error bounds when the neural-PDE (whether PINN or a Neural Operator) fails to comply with the physics. However, we believe that this is an advantage of our method. In PRE-CP formulation, the bounds are estimated across the PDE residual, where the ground truth for a well-fit solution should always be near zero. If we get wide error bars further away from the 0 for potentially high coverage estimates, it is a strong indication that statistically the model violates the physics of interest.

Consider the example with the Advection equation. We have two models, a well-fit (good model) and a poorly fit one (bad model). As shown in figure 22 in the revised paper, though we obtain guaranteed coverage in the case of both the bad and good models, the width of the error bars indicates the quality of the model. Taken for 90 % coverage, the width of the coverage bounds obtained over the bad model is substantially larger than that obtained by the good model.

There still could be a concern as to what width can be considered to be within a good range within the residual space. This could be estimated by running the PRE convolution operator(s) across a single numerical simulation of the interested physics, thereby estimating the impact of the operator in estimating the residual. The PRE over the simulation data will allow us to judge what ranges for the coverage width differentiate between a "bad" and a "good" model. Thus, PRE can be utilised as a metric to evaluate the physics deviation, and as a nonconformity score to get bounds on the physics deviation. We have added this capability to the revised paper under appendix K.