PICProp: Physics-Informed Confidence Propagation for Uncertainty Quantification
Authors: Qianli Shen
Wai Hoh Tang
Zhun Deng
Apostolos Psaros
Kenji Kawaguchi
Venue: NeurIPS 2023 poster (2023)
Keywords: physics-informed learning
uncertainty quantification
deep learning
TL;DR: We propose a method, termed Physics-Informed Confidence Propagation (PICProp) to propagate the confidence from data locations to the entire domain in terms of confidence intervals, with probabilistic guarantees.

Abstract: Standard approaches for uncertainty quantification in deep learning and physics-informed learning have persistent limitations. 
Indicatively, strong assumptions regarding the data likelihood are required, the performance highly depends on the selection of priors, and the posterior can be sampled only approximately, which leads to poor approximations because of the associated computational cost.
This paper introduces and studies confidence interval (CI) estimation for deterministic partial differential equations as a novel problem.
That is, to propagate confidence, in the form of CIs, from data locations to the entire domain with probabilistic guarantees.
We propose a method, termed Physics-Informed Confidence Propagation (PICProp), based on bi-level optimization to compute a valid CI without making heavy assumptions.
We provide a theorem regarding the validity of our method, and computational experiments, where the focus is on physics-informed learning. Code is available at https://github.com/ShenQianli/PICProp.

PICProp: Physics-Informed Confidence Propagation for Uncertainty Quantification
Qianli Shen
Wai Hoh Tang
National University of Singapore
National University of Singapore
Singapore Chenqianli@u www.iiiu.uuuueeeem edu
Singapore waihoh. waihoh.tang@nus edu. sg
Zhun Deng
Columbia University
USA
zhun. www.ddddllmaa..... edu
Apostolos Psaros
Kenji Kawaguchi
Brown University
National University of Singapore Singapore
USA a. wwwppsssbrooo..eeo edu
kenj Cnuus edu sg
Abstract
Standard approaches for uncertainty quantification in deep learning and physics- informed learning have persistent limitations. Indicatively, strong assumptions regarding the data likelihood are required, the performance highly depends on the selection of priors, and the posterior can be sampled only approximately, which leads to poor approximations because of the associated computational cost. This paper introduces and studies confidence interval (CI) estimation for deterministic partial differential equations as novel problem. That is, to propagate confidence, the form of CIs, from data locations to the entire domain with probabilistic guar- Confidence Propagation
antees. We propose a method, termed Physics-Informe (PICProp), based on bi-level optimization compute valid CI without making heavy assumptions. We provide theorem regarding the validity of our method, and computational experiments, where the focus on thysics-informee learning. Code is available at
Introduction
Combining data and physics rules using deep learning for solving differential equations (25 and learning operator mappings (20) has recently attracted increased research interest. In this context, physics-informed machine learning is transforming the computational science field in an unprece- dented manner, by solving ill-posed problems, involving noisy and multi-fidelity data well as missing functional terms, which could not be tackled before [16]. Nevertheless, uncertainty mainly due data noise and neural network wer-paraaettrization
must be considered for using such models safely in critical applications (24]. Despite extensive work on uncertainty quantification (UQ) for deep learning [1,, standard UQ has persistent limitations. Indicatively, although for most practical applications where only confidence intervals (CIs) are required, the harder problem of UQ often solved first using simplifying assumptions, and the Cls are subsequently computed. Further, strong assumptions regarding the data likelihood are required, the performance highly depends on
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

the selection of priors, and the posterior can be sampled only approximately. which leads to poor approximations because of the associated computational cost.
In this paper, we introduce and study CI estimatiin for addressing deterministic problems involving novel problem. That is, to propagate confidence, in the
partial differential equations (PDEs) form of CIs. from data locations to the entire domain with probabilistic guarantees. We name our framework Physics-Informin Confidence Propagation (PICProp) and propose practical method based on bi-level optimization to identify valid CIs in physics- thysics-informed deep learning. Our method can be applied in diverse practical scenarios, such as propagating data Cls, either estimated or given by experts, to the entire problem domain without distributional assumptions. Indicative applications are provided in Appendix [A]
PICProp provides distinct advantages over traditional UQ approaches in the context of physics- informed deep learning. Specifically, PICProp does not rely on strong assumptions regarding the data likelihood and constructs theoretically valid Cls. Further, PICProp offers a unique capability construct joint CIs for continuous spatio-temporal domains. In contrast, UQ methods typically construct marginal intervals for discrete locations. Furthermore, due to the conservative propagation of uncertainty, PICProp produces conservative CIs, which is preferred in critical and risk-sensitive applications 7 33,, considering the significant harm of non-valid confidence intervals. This sets apart from commonly used UQ methods like deep ensembles [8 and Monte-Carlo dropout [8 10
which are known to be overconfident in practice (7.
The paper is organized as follows: In Section (and its Appendices B Ch. we introduce the PDE problem setup and discuss related work. In Section 3 (and its Appendices F). we introduce the problem of CI estimation for PDEs and extensively describe joint solution predictions. Next, formulate this novel problem as a bi-level optimization problem in Section[aand Appendices G and adopt various algorithms for computing the required implicit gradients Appendix Finally, evaluate our method on various forward deterministic PDE problems in Section 5 and Appendices
and discuss our findings in Section 6
2 Background
2.1 Problem Setup
Consider the following general PDE describing a physical or biological system:
Flu(2), X(x)] = 0, x E Ω
Blu(x), [(x)] = b(x),
(1)
where x is the D dimenssiia space-time coordinate, Q is a bounded domain with boundary I, and F well B are general differential operators acting Ω and r, respectively. The PDE coefficients are denoted by (xx), the boundary term by b(x), and the solution by u(x).
In this paper, we focus the data-driven forward deterministic PDE problem, where (i) operators F, and coefficients À are deterministic and known; (ii) datasets D, { for the force condition and Dr {(=r,, ر for the boundary condition are given; (iii) uncertainty only occurs in the boundary term due to data noise, i.e., b² = b(xz) + e for = 1,...,., NN. and the rest of the dataset is considered to be deterministic without loss of generality. We denote the function obtained by utilizing the dataset z {D,, Dr} for approximating the solution u(x) of Equation with data-driven solver A as (x;
Clearly, different datasets 2 contaminated by random noise lead to different solutions u(2; 2) and thus, to distribution of solutions each It follows that the solution u(z; Z) can be construed as stochastic, where is used to denote the dataset as a random variable and to denote an instance of Z. In this regard, this paper be viewed as effort towards a universal framework, with minimal assumptions, for constructing CIs for solutions of nonlinear PDEs.
2.2 Related Work
Physics-Ifformed Neural Networks The physics-informee neural network (PINN) method, devel- oped by (25), addresses the forward deterministic PDE problem of Equation (T) by constructing a neural network (NN) approximator up(x), parameterized by 0, for modeling the approximate solution

---

u(2; The approximator (x) is substituted into Equation 1 via automatic differentiation for producing fa F Flug A] and by Blug, Al. Finally, the optimization problem for obtaining in the forward problem scenario given a dataset is cast as
0 0= argmin Ppinn (0,z), where Lpinn
ff N
Wb No
(1bo(zz))
by 112 (2)
and {wy, } are objective function weights for balancing the various terms in Equation 2 see 23
Following the determination of 0, ug can be evaluated at any in Ω. PINN can be scalled up for high-dimensional PDEs [1].
Uncertainty Quantification for Physics-Informed Learning Although UQ for physics-informed learning is not extensively studied, the field gaining research attention in recent years. Indicatively,
37 employed dropout to quantify the uncertainty of NNs in approximating the modal functions of
stochastic differential equations; [36] presented a model based on generative adversarial networks
(GANs) for quantifying and propagating uncertainty in systems governed by PDEs; 35 developed a new class of physics-informed GANs to solve forward, inverse, and mixed stochastic problems: [ further studied Wasserstein GANs for UQ in solutions of PDEs; 34] proposed Bayesian PINNs for solving forward and inverse problems described by PDEs with noisy data. Although this field is always evolving, we refer to it in this paper as "standard" UQ in order to distinguish it from our proposed method in Section In Appendices
we summarize the limitations of standard UQ and we direct interested readers to the review paper 24 for more information.
C
Interval Predictor Model (IPM) The Interval Predictor Model (IPM) in regression analysis has a similar motivation our approach. A function of interest is approximated by an IPM and the bounds on the function are obtained. (26) proposed a technique using a double-loop Monte Carlo algorithm
to propagate a CI using the IPM. The technique requires sampling epistemic and aleatory parameters from uniform distributions to estimate the upper and lower bounds of the IPM. Our approach is different for two reasons. First, we fit PINN model directly and thus, sampling epistemic parameters may not be feasible. Second, we do not require a sampling procedure. Instead, we cast the CI estimation as bi-level optimization problem and solve it using hyper-gradich methods.
3 Confidence Intervals for Solutions of Nonlinear Partial Differential Equations
In this section, we provide the definition of CI for the exact PDE solution and introduce our main assumption.
Definition 3.1 (Joint CI for PDE Solution). A pair of functions (Lz, Uz) R x R is a p-CI for the exact PDE solution uxx) if the random event & := ((Lz(x) & u(z) S Uz(2), Vr E ft holds with probability at least p.
Note that this interval refers to the exact solution u(z). Therefore, a novel problem is posed in Definition 3.1 pertaining to propagating randomness, in the form of Cls, from the data locations, typically on the boundary, to the rest of the domain. The subscript in (Lz, U2) signifies that the CI is random as it depends denoted as 2, or a collection of instances; see Appendix [ for the considered problem scenarios. Such instances can, for example, be drawn from repeating independent experiments. The subscript Z in (Lz, Uz) is dropped if there is no ambiguity.
the random dataset Z. The available data can be instance of Z,
We also note that the joint CI (also referred to as confidence sequence 17 32 in Definition 3.
corresponds to joint predictions u(-) over the whole domain Ω, rather than to marginal predictions u(x) for given The importance of joint predictions emphasized in recent works for active learning [30], sequential decision making [33], and supervised learning 222] Nevertheless, the present paper is the first work to consider joint predictions, and general CI estimation, in the context of physics-infomee informedd learning. In Appendix B we provide several theoretical remarks delineating the differences between marginal and joint predictions, as well as unions of bounds, which are typically used in standard UQ In Appendix [F] we highlight these differences with a toy example.
essential component in UQ has been
Next, to relate the PINN solution with the exact PDE solution u(x) of Definition 3.1 following assumption regarding the noise at the boundary.
we make the
3

---

Assumption& 3.2 (Zero-mean Noise). The noise on the boundary is zero-mean, i.e., Ble²] = 0, Vj ......,N):
Next, clean dataset (without noise), denoted E, can be considered the expectation of Z, i.e.,
[ In this regard, consider a data-driven PDE solver A, which yields an approximate solution u,((2; given dataset 2. We define the properness of A based on the concept of a clean dataset. Definition 3.3 (Proper PDE Solver). A PDE solver A is 1) 7-proper with respect to a clean dataset ž, if max,en [2; 2) u(x)] S 17.
Henceforth, we only consider proper PDE solvers. That is, powerful solver is available such that the absolute error of the solution obtained by the solver with a clean dataset is bounded. In all cases, assume that such an informative clean dataset exists, but is not available. Further, the constant 7). which is only related to E and A, is considered known and fixed in the computational experiments. Note that assumptions are made about u(x; 2) for * 2, i.c. there are assumptions about the PINN solver robustness against noisy data. Finally, in all NN training sessions we consider no NN parameter uncertainty, although it is straightforward to account for it by combining method with Deep Ensembles [18), for instance. The aim this paper is to propose a principled and theoretically guaranteed solution approach and thus, we refrain from applying such ad hoc solutions in the computational experiments and plan to consider parameter uncertainty rigorously in a future work. Nevertheless, we provide in Appendix [ comparison with Bayesian PINNs for completeness.
4 Method
In this section, we introduce our PICProp method and provide the accompanying theorem, as well as practical implementation algorithms. Further, we discuss the computational cost of our method and propose modification for increasing efficiency. Specifically, based on Definition 3.1 we address the problem of propagating the uncertainty from the boundary arbitrary locations domain. Equivalently, this relates to obtaining the corresponding p-CIs, i.e., obtaining the functions (L, U) E R² x R in Definition 3.1
f2 of the
4.1 Physics-Informe Confidence Propagation (PICProp)
Our method comprises two steps: (i) p-CI 2p for the clean data constructed. (ii) Σ, is propagated to any point x E Ω of the domain by solving the following problem,
L(x) min uA(2;)) U(z) = max (A(=;)) + =EZ,
zE
which searches for configurations of points within the interval E, such that the predicted values for u(x:) given by the 77-proper solver A are minimized or maximized, respectively.
The CI in the first step can either given or constructed from data. The former case relates to various practical scenarios (e.g., the error ranges of measurement sensors with certain confidence degree are known or given by experts), whereas the latter serves as a data-driven scheme, where construct the CI in the first step depend on the data
less prior knowledge is required. Techniques (noise) distribution. For example, if the noise is assumed to be Gaussian, the corresponding CI can be identified by chi-squared Hotelling's 's T-squared statistics. If the noise distribution is not known, concentration inequality can be used to identify the CI. Relevant details are provided in Appendix [ See also Figure for illustration.
The interval constructed via Equation (3) is a joint p-CI for the whole domain of u. This fact is formalized in the next theorem (see Appendix for the corresponding proof).
Theorem 4.1. Consider a p-C1 for 2 denoted as Zp.. constructed such that Pr(z E 2,) 2 p. Then the functions UL obtained by solving Equation (3) define a p-CIffo for u.
Although the robustness of the PDE solver does not affect the validity of Theorem 4.1 1. it nevertheless
affects the quality of intervals constructed via Equation (3). Specifically, although valid CIs can always be constructed with arbitrary proper PDE solvers by solving Equation (3). a robust PDE
solver less sensitive to data noise is expected to yield tighter CIs compared to a less robust solver.

---

Note that the confidence guarantee of the constructed CIs corresponds to probability that is greater than, rather than equal to, p. which means that the constructed CIs are conservative. From design perspective, conservative uncertainty estimation is the preferred route for critical and risk-sensitive applications 13,, although it can more costly.
Next, a brute force and computationally expensive approach to solve Equation (3) pertains to traversing E žp in an exhaustive search (ES) manner and invoking PDE solver for each candidate z. Alternatively, we propose first employ a NN for obtaining 4d(x) by minimizing the physics-
informed loss of Equation (2 and approximating u u(2; 2) for each Zp. In this regard, if PINN is assumed to be an 17-proper PDE solver, Equation can be further extended to a bi-level optimization problem cast
L(x)= = min U8(x) (x) - 77 U() = max U8(z) (I) + 17. zEZp
=E2,
s.t. 8(z) arg min Cpinn (0, 2).
Nevertheless, explicit gradient V. calculation is intractable, multiple steps of gradient descent for the NN are required for a reliable approximation. Despite this fact, bi-level optimization involving large-scale parameterized models such NNs, with main applications in meta-learning, has been extensively studied in recent years and has resulted to numerous implicit gradient approximation methods (see Appendix D In this regard, Algorithm in Appendix H summarizes our method for solving efficiently Equation (3) by using such approximations. In passing, note that r is considered for all implementations this work.
4.2 Efficient Physics-Informne Confidence Propagation (EftiPICProp)
If joint Cls are required for multiple query points, denoted by the set Xq. Algorithm T has to be used multiple times to solve Equation (4) for each and every point the set. This clearly ineurs a high computational cost. To address this issue, straightforward approach pertains to training an auxiliary model parameterized by & that models the as (I,) (x,)) for all query points T. E Xg or more generally of the entire domain Q. However, this naive regression using solutions of the bi-level optimization at a set of locations, which we refer to as SimPICProp, does not provide accurate results when the size of the query set is small; see computational results in Section 5.1
For improved generalization properties with limited query points, we propose train a two-input parameterized meta- model ug (Ig. .) that models the resulting PDE solution u(x), via Equation (4) for each different query point. For example, note that T.(x:) Equation (4) evaluates the NN (ā(=) at the location of the query point and z minimizes the corresponding value. As a result, a trained model (xg. minimizes the value of at Zg. the sought y-CI can be simply obtained by evaluating
uy at (Ig. Tg). More generally, given a query set consisting of K randomly selected query points
X, = [ the model is trained such that it solves the problem
min
Σ Cpinn (1 : ,), (xt)) + Lpinn (uy (a 2), (22(1)1))
2K
where 2 (a arg mig Uā(a) (x))), ),, 2 (a = arg max u Uz(z) (1 zE2,
(5)
zE2,
s.t. ((=) arg min Lpinn (ua (-)())z2)
The estimated p-CI is, subsequently, given by
L(x) - u((, U(x)= = u (x)
(6)
In fact, we can train a single model for both upper- and lower-bound predictions because the parameters shared in Equation (5). This is achieved by encoding an indicator into the input vector, so that u (Zg. = up (Zg. - II)and and u (Ig.:) us ' 1). Furthermore, to balance the weights of the terms in the objective Equation (5) more effectively, PINNs are separately trained with boundary conditions z (z and 2 22 In this regard, Uy is obtained by
5

---

K
min Cem (ig (x):K;0 ; where
Lem ==
1 2K
() l 15) +)+
E re
( 27 2)
s.l. L(syy 1;) = (up( -1) "0(=* (x)) + (ubb (z,, x. "ă (xf)) (x:))
and A [0 is a weight to balance the two loss terms, of which the first aims enforce close predictions CIs at query points, and the second to utilize extra information in PDE solutions for better generalization. We refer to this model as EfriPICProp and provide the corresponding implementation in Algorithm 2 in Appendix H Specifically, the EffiPICProp version with = 0 is called SimPICProp. as it simply focuses on close predictions of confidence intervals at query points, while disregarding additional information at other locations.
5 Computational experiments
We provide three computational experiments that are relevant to physics-informet machine learning. We will start with ID pedagogical example illustrate our proposed approaches. Next, we will move to more complicated problems involving standard representative equations met in many real-world problems including heat transfer, fluid dynamics, flows through porous media, and even epidemiology 21 (88. Details regarding implementation, query points, and computational times can
be found in Appendix [
5.1 A Pedagogical Example
We consider the following equation, which has been adopted from the pedagogical example of (1].
Ugr uuu U2 f(x),
[-1,1], (PDE),
where fxx = -x² sin (Tæ) - π cos((TT) sin (IT.). For the standard BCs u( = u(1) 0, the exact solution is u(x) sin( The purpose of this one-dimensional PDE with only two boundary points, i.e., with a two-dimensional search space in Equation (4), is to demonstrate the applicability of method in different data scenarios (see Appendix [D. Further, the small size of the problem is ideal for comparing computationally expensive methods, such as exhaustive search (ES) and employing PICProp separately at all domain locations, with the more efficient implementation of Algorithm 2
(a) x2 CI
(b) CI
(c) Hoeffding CI
0.2
x² CI sample
Cl samples
Hoeffding
0.1
ofx= 11
-0.1
0.2
-0.1
0.0
0.1
0.2
0.2
-0.1
0.0
0.1
0.2
0.2
0.1
0.0
0.1
0.2
uix 1)
uix
1)
u(x
11
Figure 1: Pedagogical example: Boundary 95% Cls to be propagated to the rest of the domain. (a-b) Random samples follow N(0, I²I), where a = 0.05. The black crosses correspond to 5,000 random samples from the distribution. A x2 CI is constructed using the only one available sample and T CI is constructed using the five available samples. (c) Random samples follow U[[-1.5, 1.5) x [-1.5, 1.5)). The black crosses correspond to 2,000 random samples from the distribution. A Hoeffding CI is constructed using the five available samples. All available samples are shown with red crosses.
We consider two types of noisy BCs (Figure D First, we consider BCs with Gaussian noise given as
u(-1) NN(0, u(1) ~ o(),
(9)

---

where the standard deviations are di 02 0.05, i.e., the data is assumed to be Gaussian and the standard deviations are known (or estimated if unknown). We sample one or five datapoint(s) for each side of the boundary, for each case, respectively. Given the dataset, x² CI and a T CI are constructed according the corresponding problem scenarios discussed in Appendix D see also Figure a)) ab).
Second, consider BCs with uniform noise given
u (-1 ~ Li),
~ u(1)~ D2),
(10)
where D1 = D2 1.5 are assumed to be known. We sample five datapoints for each boundary side. Given the dataset, a Hoeffding CI constructed; see also Figure T.).
For the query set, six points are used for implementing SimPICProp and EffiPICProp in Section 4.2 whereas 41 points used for implementing the brute-force version of PICProp, which amounts to solving the bi-level optimization problem of Equation (4) for each point. The implementation details of the PICProp-base based methods are summarized in Table J.1 Further, 5,000 BCs are sampled for solving Equation 4 with ES, which can be considered as reference solution.
Chi-squared
T-squared
Hoeffding
1.0
ETRPICProp
EMPICProp PICPhop
RCPra
PtEPop
0.5
ES
ES sample
ES
sarT ple
sample
0.0
-0.5
-1.0
-0.5
0.0
1.0
-0.5
1.0 00.5
0.0
0.5
1.0
Figure 2: Pedagogical example: CI predictions for the entire domain r [-1, Left: Noise is Gaussian with known standard deviation, and x2 CI is constructed the boundary. Middle: Noise is Gaussian with unknown standard deviation, and T CI is constructed at the boundary. Right: Uniform noise with known bounds, and a Hoeffding CI constructed the boundary.
The CI predictions for the three cases are visualized in Figure 2 whereas the meta-learning curves be found Appendix K.1 Note that the sizes of the CIs are approximately equal to 0.2 near the boundaries, which are comparable with the propagated CI sizes of Figure [ and become larger around æ of +0.5.. In all cases, the PICProp results match those of ES, despite the fact that PICProp requires solving a challenging bi-level optimization method with approximate implicit gradient methods. Further, the efficient method EfiiPICProp performs accurately as compared its reference method PICProp, while retaining significantly lower computational cost (see Figure for x² results and Figure K.5 for 2 and Hoeffding results). On the other hand, the predictions of SimPICProp match the brute-force PICProp results only on the query points, but generalize poorly on other locations. contrast, EffPICPrrp, which is based on meta-model that fits the PDE solutions for multiple query points, generalizes well unseen locations.
Some additional experimental results summarized in Appendix K.1 for the completeness of empirical study: Firstly, a comparison among PICProp and EffiPICProp with various As is provided in Figure K.3 to better demonstrate how the balance between two EffPICProp loss terms in Equation 7 affects the CI prediction. Next, we conduct an additional experiment to empirically verify the validity of EffPIICProp CIs Appendix K.1.4 Finally, brief comparison with Bayesian PINNs presented in Appendix
5.2 Two-dimensional Poisson equation
Next, consider two-dimensional Poisson equation given as
Ugz Uhy f(x (I. e u(xx 3) = e" te I = (BCs).
(PDE),
(11)

---

Chi-squared
Chi-squared (zoomed in)
1.2
1.0
REFrp
semple
0.5
Parta EIRRChop
0.8
0.0
0.6
-0.5
Serfichep
0.4
-1.0
Farb STAROPO
-1.0
0.5
0.3
0.6
Figure 3: Pedagogical example, x²: Although the results of both SimPICProp and EftiPICProp match PICProp results on the query points, EfriPIPProp generalizes better on unseen locations.
The exact solution u(z, = e", and the corresponding f(z y) is obtained by substituting
u(x, y) into Equation [11). This example serves to demonstrate the effectiveness of our method for treating cases with large numbers of noisy BC datapoints. such cases, an exhaustive search approach for solving the bi-level optimization of Algorithm generally fails to identify the global optimum, unless we increase the number of trials and the computational cost at prohibitive levels.
We propagate the following fixed CI from the boundary to the rest of the domain (ær, v) E (-1,1) x 1)1
u(x,y) E [e" 6" + + + 0.05), I +r
(12)
To solve the problem, 10 boundary points are sampled for each side of the domain, i.e., 40 boundary points total are used for training and uncertainty is considered at all points. The search space is thus 40-dimensional. The implementation details of the PICProp-based methods are summarized in Table J.T whereas 4,000 steps of standard PINN training fair comparison, 121 query points for PICProp are selected in the domain, see Figure [.1 and the computational time is approximately the same as of ES with 1,000 trials. Moreover, include the ES result with 5,000 trials, which is considered to be exhaustive, as reference.
for each BC for ES. For
(a.1) = 0.5
2) 0.5 (zoomed in)
(b.1) 0.0
(b.2) 0.0 (zoomed in)
3.0
3.9
3.5
2.3
3.0
2.2
3.8
2.0
1.3
3.5
2.0
1.9
1.0
1.5
-0.5
0.0
3.3
-0.
-0.1
-1.0
-0.
0.0
1.7
-0.2
-0.1
(zoomed
(zoomed
3.0
4.1
5.5
2.0
6.0
3.9
3.5
2.7
3.8
2.6
2.5
30
2.4
35
2.0
3.0
-1.0
-0.3
-1.0
Exact
EHtiPICProp (121 points)
ES (1000 points)
ES (5000 points)
Figure 4: 2-D Poisson equation: Slices of CI predictions for (-0.5, 0.0, 0.5, can identify the solution to the bi-level optimization problem more efficiently and accurately, as compared with the computationally expensive, search-based method ES.
0.5, 1.0). EMpPICPrrp
Figure plots the slices of CI predictions includes slices along the y axis and the heatmaps over the entire domain. For reference, the size
E {-0.5, 0-0, 0.5, 1.0], whereas Appendix K

---

of the CIs obtained with EffiPICProp is approximately 0.1 across the domain, which is equal to the size of the propagated CI of Equation (12). The predicted CI size is also larger in most cases than that of the ES-based Cls. We conclude that 1. 1,000 trials for ES are clearly not sufficient, as the intervals of ES with 5,000 trials are significantly wider; thus, a global optimum is not reached with 1,000 trials. In contrast, EffiPIPProp converges to similar or wider intervals compared to ES with 5,000 trials. Overall, the optimization-ben baaed method Prri can identify the solution to the bi-level optimization problem more efficiently and accurately, as compared with the computationally expensive, search-based method ES.
5.3 One-dimensional, time-dependent Burgers equation
In this example, we study the Burgers equation given as
(PDE), u( t 0) = (TT), u(x= = -1,)) u(z 1,t) 0, (BCs),
Ue + uz Vuga 0, (x. t) =|-1, [0,
(13)
where = 0.01 is the viscosity parameter. We consider clean data for the BCs and noisy data for the i initial conditions (ICs). Specifically, for the ICs, we propagate the following fixed CI from = 0 to the rest of the time domain E (0, 1] and for all [-1,
u(x f = 0) E [- -in((TT) 0.2, - sin(Txx + 0.2].
(14)
This problem serves to demonstrate the reliability of our method even in challenging cases with steep gradients or discontinuities of the PDE solution. For solving the problem using PINNs, we consider 200 boundary points, 256 initial points, and 10,000 samples in the interior domain. The 256 initial points are evenly spaced along [-1,1] 1]. The implementation details of the PICProp-based methods summarized in Table whereas 2,500 and 5,000 trials are considered for ES. The 246 PICProp query points are visualized in Figure J.I and the respective computational time is approximately equal that of ES with 2,500 trials.
Slices of CI predictions for E {0.0, 0.25, 0.5, 0.75} are plotted in Figure whereas Appendix K.3 includes the respective meta-learning curves. We conclude that EffiICPrrp, although computa- tionally cheaper, in most cases provides wider, i.e., more conservative, Cls, as compared to ES with varying numbers of trials. It also succeeds in capturing the expected increase in uncertainty as the query points move away from the source of uncertainty, which corresponds to larger values. Specifically, the size of the CIs propagated CI in Equation (14), and becomes approximately 2 at ( = 0.75 and near x = 0. On the other hand, since we consider 256 random initial points, the search space of ES is 256-dimensional. It requires a prohibitively large number of trials to identify the global optima of Equation 4 Finally,
0 is approximately 0.4, which is equal to the size of the
note that the results of EMiiICPPrp in Figure ] to 0 in Equation 7 because we found that the additional information provided by the second loss term in Section 4.2 adversely affects generalization. This is partly due to the large number of query points in this problem.
6 Conclusion
Physics-informet machine learning is transforming the computational science field in an unprece- dented manner, by solving ill-posed problems, involving noisy and multi-fidelity data as well as missing functional terms, which could not be tackled before. However, uncertainty quantification necessary for deploying such methods in critical applications. In this work, we considered a novel problem pertaining to propagating randomness, in the form of Cls, from data locations (at the boundary but can easily be extended to inner or collocation data points) to the entire domain, with probabilistic guarantees. We proposed PICProp. method based on bi-level optimization, and its efficient version, with a theorem guaranteeing the validity of the obtained CIs (Algorithms I and Theorem 4.1. To the best of our knowledge, this is the first work joint CI estimation with theoretical guarantees in the context of physics-ifformed learning, while our method can be applied in diverse practical scenarios without imposing any distributional assumptions.
In the computational experiments of Section S involving various partial differential equations, we demonstrated the effectiveness and efficiency of our proposed methods. Specifically, the predicted CIs are wider, i.e. more conservative, in most cases, as compared with the more computationally demanding exhaustive search for solving the bi-level optimization. Although not encountered in the

---

(a.1) 0.0
(a.2) Izopmed in)
0.25
(b.2) 0.25 (zopmed
044
0.5
0.5
0.5
5-5
-05
-05
-10
-10
0.5
(c.2)
(d.1)
1.0
10-
1.0
-0.5
-0.5
-10
Exact
EMiPICProp (246 points)
ES (2500 points)
E5 (5000 points)
Figure 5: Burgers equation: Slices of CI predictions for {0.0, 0.25, 0.5, 0.75). EfiPICProp. although computationally cheaper. in most cases provides wider, i.c.. more conservative Cls, as compared to ES with varying numbers of trials. interpret the cost ES, note that its search space is 256-dimensional and requires prohibitively large number of trials to identify the global optima of Equation
experiments, potential limitation of our method is that the obtained CIs can be loose in some cases, even when propagating tight from the boundaries, due to the max/min operations in the bi-level optimization of Equation (3). Because of the theoretical challenges involved, we plan to study this limitation and perform pertinent experiments a future work.
Acknowledgements
This material is based upon work supported by the Google Cloud Research Credit program with the The computational work for this article was partially performed award resources of the National Supercomputing Centre, Singapore /htpss//www.nsccc.gg)
References
[1] Moloud Abdar, Farhad Pourpanah, Sadiq Hussain, Dana Rezazadegan, Li Liu, Mohammad Ghavamzadeh, Paul Fieguth, Xiaochun Cao, Abbas Khosravi, U Rajendra Acharya, et al. A review of uncertainty guantification in deep learning: Techniques, applications and challenges. Information Fusion, 76:243-297, 2021.
[2] Vidmantas Bentkus. An extension of the hoeffding incquality to unbounded random variables. Lithuanion Mathematical Journal, 48(2):137-155) 2008
[3] BVSS Bharadwaja, Mohammad Amin Nabian, Bharatkumar Sharma, Sanjay Choudhry, and Alankar Alankar. Physics-informed machine learning and uncertainty quantification for mechanics of heterooeeneous materials. Integrating Materials and Manufacturing Innovation, pages 1-21, 2022.
Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng, editors. Handbook of Markov Chain Monte Carlo. Chapman and HalVCRC, may 2011.
[5] Nicolas Brosse, Alain Durmus, and Eric Moulines.
The promises and pitfalls of stochastic gradient langevin dynamics. Advances in Neural Information Processing Systems, 31, 2018.
[6] Min Cai, George Em Karniadakis, and Changpin Li. Fractional seir model and data-driven predictions of covid-19 dynamics of omicron variant. Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(7):071 101, 2022.
[7] Kamil Ciosek, Vincent Fortuin, Ryota Tomioka, Katja Hofmann, and Richard Turner. Conservative uncertainty estimation by fitting prior networks, In International Conference on Learning Representations. 2020.
10

---

[8] Morris L Eaton and Michael D Perlman. Concentration inequalities for multivariate distributions: I. multivariate normal distributions. Statisticx & probability letters, 12(6):487-504, 1991.
[9] Andrew YK Foong, David R Burt, Yingzhen Li, and Richard Turner. Pathologies of factorised gaussian and me dropout posteriors in bayesian neural networks. stat, 1050:2, 2019.
[10] Yarin Gal and Zoubin Ghahramani. Dropout a bayesian approximation: Representing model uncertainty in deep learning. international conference on machine learning, pages 1050-1059, PMLR, 2016.
[11] Yihang Gao and Michael K Ng. Wasserstein generative adversarial uncertainty quantification in physics- informed neural networks, arXiv preprint arXiv:2108.1300 2021.
[12] Riccardo Grazzi. Luca Franceschi, Massimiliano Pontil, and Saverio Salzo. On the iteration complexity of hypergradient computation. In International Conference on Machine Leurning, pages 3748-3758, PMLR. 2020.
[13] Bobby He, Balaji Lakshminarayanan, and Yee Whye Teh. Bayesian deep ensembles via the neural tangent kernel. Advances in neural information processing systems, 33:1010-1022, 2020.
[14] Zheyuan Hu, Khemraj Shukla, George Em Kamiaaakis, and Kenji Kawaguchi. Tackling the curse of dimensionality with physics-informet neural networks. arXiv preprint arXiv:2307. /2306, 2023.
[15] Kaiyi Ji, Junjie Yang. and Yingbin Liang. Bilevel optimization: Convergence analysis and enhanced design. In International Conference on Machine Learning, pages 4882 4892. PMLR, 2021.
[16] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris. Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422 2021.
[17] Tze Leung Lai. On confidence sequences. The Annals of Statistics, pages 265-280, 1976.
[18] Balaji L.akshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in neural information processing systems, 30, 2017.
[19] Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing millions of hyperparameters by implicit diffcrentiation. In International Conference on Artificial Intelligence and Statistics, pages 1540-1552.
PMLR, 2020.
[20] Lu Lu. Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identi- rying differential equations based on the universal approximation theorem of operators. arXiv preprint 1xXiv:110.031193, 2019.
[21] Xuhui Meng, Hessam Babaee, and George Em Karniadakis. Multi-fidelity bayesian neural networks: Algorithms and applications. Journal of Computational Physics, 438: 10361, 2021.
[22] Ian Oshand, Zheng Wen, Seyed Mohammad Asghari, Vikranth Dwaracherla, Botao Hao, Morteza Ibrahimi, Dieterich Lawson, Xiuyuan Lu, Brendan O' O'Donoghue, and Benjamin Van Roy. The neural testbed: Evaluating joint predictions. arXiv preprint arXiv:21 1n..0629, 2022.
[23] Apostolos Psaros, Kenji Kawaguchi, and George Em Karniadakis. Meta-leaming pinn loss functions. Journal of Computational Physics, 458:1 1 121, 2022.
[24] Apostolos F Psaros, Xuhui Meng, Zongren Zou, Ling Guo, and George Em Karniadakis. Uncertainty quantification in scientific machine lcarning: Methods, metrics, and comparisons. Journal of Computational Physics, page 111902, 2023.
[25] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
[26] Jonathan Sadeghi, Marco De Angelis, and Edoardo Patelli. Robust propagation of probability boxes by interval predictor models. Structural Safety, 82:101889 2020.
[27] John G Saw, Mark CK Yang. and Tse Chin Mo. Chebyshev inequality with estimated mean and variance. The American Statistician, 38(2):130- 132, 1984.
[28] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, llya Sutskever, and Ruslan Salakhutdinov. Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research,
151): 2014.
[29] Bartolomeo Stellato, Bart PG Van Parys, and Paul J Goulart. Multivariate chebyshev inequality with estimated mean and variance. The American Statistician, 71(2):123-127, 2017.
[30] Chaoqi Wang, Shengyang Sun, and Roger Grosse. Beyond marginal uncertainty: How accurately can bayesian regression models estimate posterior predictive correlations? In International Conference on Artificial Intelligence and Statistics, pages 2476-2484, PMLR, 2021.
[31] Zhuo Wang. Chen Jiang. Pengwei Liu, Wenhua Yang. Ying Zhao, Mark F Horstemeyer. Long -ingQQig Chen, Zhen Hu, and Lei Chen. Uncertainty quantification and reduction in metal additive manufacturing. Npj Computational Materials, 6(1):175, 2020.

---

[32] Ian Waudby-Smite and Aaditya Ramdas. Confidence sequences for sampling without replacement. Advances in Neural Information Processing Systems, 33:20204-2022 2020.
[33] Zheng Wen, lan Osband, Chao Qin, Xiuyuan Lu, Morteza Ibrahimi, Vikranth Dwaracherla, Mohammad Asghari, and Benjamin Van Roy. From predictions to decisions: The importance of joint predictive distributions. arXiv e-prints, pages arXiv-2100 2021.
[34] Liu Yang. Xuhui Meng, and George Em Karniadakis. B-pians: Bayesian physics- informed neural networks for forward and inverse pde problems with noisy data. Journal of Computational Physics, 425: 109913,,, 2021.
[35] Liu Yang. Dongkun Zhang, and George Em Kariaaaakis Physics-informed generative adversarial networks for stochastic differential equations. SIAM Jouraal on Scientific Computing, 42(1)) 2020.
[36] Yiho Yang and Paris Perdikaris. Adversarial uncertainty quantification in physics-informed neural networks. Journal of Computational Physics, 394:136-152, 2019.
[37] Dongkun Zhang, Lu Lu, Ling Guo, and George Em Karniadakis. Quantifying total uncertainty in physics- informed neural networks for solving forward and inverse stochastic problems. Journal of Computational Physics, 397:108850, 2019.
[38] Zongren Zou, Xuhui Meng, Apostolos F Psaros, and George Em Karniadakis. Neuraluq: A comprehen- sive library for uncertainty quantification in neural differential equations and operators. arXiv preprint arXiv:2208. 2022.
12

---

A Indicative applications
Our proposed method, delineated in Algorithms and of of Appendix provides way to quantify how the data uncertainty affects the obtained solution. In addition, the unseen clean data lies within the constructed/given data bounds, the ground truth solution lics within our predicted Cls. Some indicative applications include:
In additive manufacturing problems, understanding, through our computed Cls, how the model responds to different material properties arising from the manufacturing process can help in model correction and thus quality control (see (3I).
2. Further, in material science, cracks affect the elasticity constants which together with geometric properties affect deformations (sec B). The range of these deformations can be identified with our method. Accurately predicting the upper bound of such deformations can help designers adjust cross-scctional properties structural elements to better account for potential cracks.
More generally, our model can extended to account for PDE model uncertainty. For example, PINNs have been used for short-term predictions of infected cases and deaths due to COVID-19 (see (6). In such problems, our approach can be used to provide bounds of uncertainty for the predicted variables of interest. Accurately predicting the upper bound of infected cases help the states better prepare for short-term hospitalization needs.
B Additional details regarding standard UQ
Standard UQ considers that the data is generated from a distribution p(u/ar, 0 = 8*), where 0* is the NN parameters that assumed to represent exactly the PDE solution u(x)) Given dataset z, standard UQ approximates the solution ug- (x) with the distribution
p(ula:; 2) = p(ulz, (p(00))00
(15)
where p(0|2) xc p(0)p(=(0) according to Bayes' thcorem and p(z[0) is the data likelihood. However, as a closed-form p(0)z) typically intractable, UQ draws approximate samples from p(0)z), denoted as 1 Next, it constructs p(u|a:; =) for each based on the collected samples to approximate p(u lal;; z). Typically, a Gaussian distribution is fitted to the PINN solutions {ug, [z) No and its corresponding statistical intervals are used in lieu CIs for the PDE solution. Loosely speaking, statistical intervals of (ي,} can be viewed as credible intervals from Bayesian statistics, and intervals of {Ma, (m) as the corresponding prediction intervals.
Alternatively, one could use a concentration inequality to construct directly a prediction interval using the samples {ug, (- Nevertheless, such intervals still correspond to the Bayesian estimate (e- (D) tie rather than the exact solution u(æ), unless further assumptions are made. Even if the assumption that the Bayes estimate corresponds to 8 is made, i.c., 8" = pf, the typically considered ( as ug can be arbitrary nonlincar (see also Appendix C)
is not guaranteed to equal
Standard UQ approaches have several other limitations. For example, strong assumptions regarding the noise distribution and data likelihood are often required, and depending on the problem is notoriously difficult to tune Bayesian NNs (.) Furthermore, it is not possible to obtain/sample p( of/ 2) exactly in practical problems and thus, additional approximations are required, typically referred to as approximate posterior inference. In addition, UQ approaches for PINNs are not supported by theoretical guarantecs regarding the quality of uncertainty propagation due the physics-informed constraints (II). Finally, for the case of multiple available datasets, Appendices and the limitations Monte Carlo-based approaches.
C Additional details regarding Monte Carlo Simulation
Assume that we have m i.i.d. datasets denoted as (200p and O-proper data-driven PDE solver (soc Definition 3.3). is natural to consider solving the PDE with each dataset in a Monte Carlo manner and then using a concentration inequality to construct Cls. However, CI for the exact solution u(rr) cannot be identified without further assumptions. The reason that we only have access Elu(x; 2)], rather than u(z)) given as u(x) u(2; Z) u(æ; E(2]). considering lata-driven solutions from a o-proper solver. Clearly, the difference between the two is in the form of [ ]] # [EEX]), where f in our case is the solver. In Figure C.I we provide an illustration of this point as applied to lineat/convex/concaee functions, as straightforward consequence of Jensen's 's inequality.
13

---

(0) 100
(b) (eWee
[c) (0) sin)
1.0
1.0
1.0
Fixt #nlrxs72
0.8
0.0
0.8
/ORUNTI SP
0.8
0.6
0.5
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.2
0.4
0.6
0.00
0.2
0.4
0.6
0.0
0.4
0.6
0.8
Figure C.I: Difference between (X and f(E[X), (a) Linear function (())==,B[[(x)) =
[[E[x]): (b) Convex function fxx = (e
()) ((B(X) Concave function
f(z) = E[f(X) < f(E[)].
D Confidence Interval Construction with Data
Consider a dataset pp which contains m independent samples of an r-dimensional random vector Z. Based on the data (20)) 1z=1, we construct a CI for Z, denoted as 2p, where p is the confidence level.
Equivalently, Pr(z 3,) 2 p.
this regard, we consider the following problem scenarios:
1. The random vector Z follows a multivariate Gaussian distribution with known covariance matrix Σ of rank Given observation z, i.e., = 1. the chi-squared statistic x² = (=
1)) follows a chi-squared distribution with n degrees of freedom. The corresponding p-C1 is Σ, v(= =² 2' 2") and A:X" is the squared root of the upper (1 p)/2 critical value of the chi-squared distribution with degrees of freedom; see Figure for example corresponding the case
{='l 1/2==''|>=- < k ]. where ||= 112-=1lp ='1|=->
2. The random vector follows a multivariate Gaussian distribution with unknown covariance matrix Σ. and m observations are available. The sample mean and covariance are estimated 5 Σn and
- En (zi (2; 2)T, respectively. The Hotelling's 's t-squared statistic t [= μ1 - n) follows Hotelling's T'usquredd distribution with parameters n. m The corresponding p-Cl 2, [='1 2 o'lls,-1 S m -1/2 k...m........ where ||= 2 -1ss;-1
V (a (U'T (z 2²) and A-T k...m-1.p denotes the squared root the upper (1 - p2/2 critical value of the aforementioned Hotelling's T --squaeed distribution; see Figure for example corresponding to the case of n and m
3. The components of the random vector Z follow unknown distribution, and independent observa- tions are available. For this case, concentration inequality can used each component, i.c., for each data location, and the obtained bounds can be combined using a union bounds for obtaining a p-CI (see also Appendix E.
Note that many popular concentration incqualitics require specific assumptions for the unknown distri- bution besides the fact that independent samples are available. For example, Hoeffding's 's incquality and its variants require the random variable be bounded. [2] extends Hoeffding' inequality to unbounded random variables, but is limited Pareto-type tails. Other concentration inequalities without heavy assumptions for distribution, such as Chehyshev's inequality and Chernoff's inequality. require various true population statistics (e.g., the population mean and variance) for computing the bounds, which are not available in the cases consider. In this regard, 27 provides Chehyshev-type hound with estimated mean and variance. [29] further extends the bound of [27 the multivariate case.
4. The CI ž, can alternatively be given by domain experts.
E Joint Confidence Intervals for PDE Solutions
In this section, we provide the necessary definitions, lemmas, and remarks that are required for introducing and addressing the Cl propagation problem in Section 4.1 That is, how propagate, with probabilistic guarantees, a CI constructed at the locations of the available data arbitrary sets of query points of the domain.
To this aim, first define the marginal CI for single query point xq.
14

---

Definition E, (Marginal CT). A pair of random variables (L, E is' is p-CI for u(x;,) on the single query point Biy E 92, if the random event := (L < v(x,)) < () holds with probability least p.
Intuitively, considering repeating experiments of drawing random datasets (2(0) fios the proportion of the (L, U) Cls, constructed according to Definition E.1 that contain u(Zy) is at least p. Next, we define the joint CI for multiple query points.
Definition E.2 (Joint CI), A pair of functions (L,U) I2 x R is a p-CI for on query set X, S S2, if the random event & == Xx (L(x)) < u(xa) < (24)] holds with probability at least p.
For simplicity, in the rest of the paper we refer to the p-CI for u on the whole domain ft as p- p-Cl for u, and to the p-CI for on Xg C Q as query set p-CI.
Remark E.3. Let (I,, U) be a p-CI for For any X₂ St, (L. U) is p-CI for Xgg, i.e., the random event E(y, Z) == 2x (L(z)) < u(x,)) < (xa)) holds with probability at least p.
Remark E.3 indicates that a p-CI for on an arbitrary query set can be derived from a p-CI for u by slicing. The query set X, can be either discrete continuous including multiple even infinite queries. Note that the p-Cls of Definition
(E.zitiionnEEaaaae E.22 and Remark [E.3 Eh3ld with probability at least p; i.c., there no statement about their tightness. In this regard, the next remark shows that slicing using larger sets Xg yields tighter p-Cls for u on Xg. Nevertheless, the maximum tightness that can be achieved, for the p-CI for on Xg. is determined by the tightness of the p-CI for u on the whole domain Xg. This is corroborated by the following remark.
Remark E.4. Let (I, U) RD x Rd be p-CI for For query set X. C X2 C Q, the probabilities of random events E == Ex (L(zg) u((a) < (xxz)) &g : Nagez (L(rg) u(zr) ((zy)) satisfy
Pr(E1 2 Pr(E2).
(16)
Next, we define the union of single-query p-Cls, It is often used in practice as a proxy to the query set p- and can be used also within the context of the present study for addressing problem scenarios 3 and 4, as described in Appendix [
Definition E.5 (Union Cls). A pair of functions UL E R x R constructed based on (=(0) Jmm a union of p-CIs for on query set Mg S Ω if for each xy E Xz. (I(x,), U(1,)) is a p-CI for u(T,)), i.e., each random event & = (I(z,) < u(z,) (xx)) holds with at least probability p.
The next lemma provides insight into the connection between the marginal p-CI on X, and the union of p-Cls. Lemma E.6. For and Nz. p-CI is union of p-Cls, whereas union of p-Cls is not necessarily a p-CI.
That is, u(x,), Vag E Xy. belongs to the p-CI constructed by slicing of the p-CI for with probability at least p. On the other hand, the probability of the event Oz, Exx, (L(ag) u(x4) < Uxx)).. where (L, U) are constructed based well as the computational experiment of Appendix
union p-CIs, may be less than p. This supported by the following two remarks as
Remark E.7. The probability of the random event Ex, (L(,)) u(2,) U(2,)), }. where (L. are constructed based on union of p-CIs, is upper-bounded by the minimum probability of the events (L < u(r,)) < (24)), Va, E Xg. That is,
Pr(Ni,ex, (M(,)) < u(za) < ((x)))
S min Pr( (1))) < <u(s)) < U(z)))). a,EX;
(17)
Remark E.S. If one of the p-CIs is tight, i.e. Enly E X, such that Pr((L(u)) < u(x)) < ((())) union of these p-CIs is a p-CI for if and only if holds Va,, E X2 (m%) that U(x1)))= 1.
p. the
Pr([L(x,,) < u(x,)) < ((= < u(z)))
(18)
Finally, in order obtain guarantee that union of CIs is p-CI, we can equate lower bound of the probability the former with p. This lower bound can obtained by disregarding correlations among the different locations, as formalized in the following lemma.
Lemma E.9. For u and finite query set X,, a union of (1 - (1 -11 (1-p))/x))--Cssss p- u as
Pr(nz,ex, (L(x,) u(x,,<()xx)))) =1 - PM(uxxx (Lxx) < f(x) 0)
21 - Σ (1 (L(( A u(s)= s ((x))))
(19)
1- P ) = p- [x,
21 - 1X,111 - (1
15

---

However, as the size of the query set X, increases, the required confidence level of each interval must approach which leads loose union bounds. Finally, another approach for obtaining p-CI for multiple queries is use multivariate concentration inequalities []] Nevertheless, this approach limited a finile number of queries.
F Toy Linear System
In this section, we highlight the differences between joint CIs and union of Cls, described in Appendix H
with simple example. Consider function and model given uo()) (1 - 2:)0 + 201, [0. 1T. with (fo: ) ?R?, We assume that there exists such that ug-- = u, and that we are given measurements
(Zo, Z1) (0 1). Thus, one observation Z. we simply estimate i Z. As result, the estimator
for is given as
0g x) Zo =21 N =)0 201, 2x
2x 1).
(20)
Next, define intettrval the lower and upper intervals (Le, Ux) obtained by subtracting and adding k standard deviations from the estimated of u, respectively, i.c.,
(x(x) == - + x71 KV 2x
(21)
(x(z)= = =0 + x/1 + kv 2x 2x 11.
Remark F.I. We define Ky -interval as the interval constructed a union of one-query p-Cls for uo uo(x) and for different values of r Q. similarly to Definition E.5 The value kp the upper (1 - p2/2 critical value of the standard normal distribution.
Proef: For x E Ω
Pr(1Tay S uxxx)= s Uny})
=Pr(1(1 - x) xbi - Kn /2x² 2x +1 f(x) < x) + xb1 + k2 V222 2x +1+1 (22)
-Pr([|u(T)) - (1], < K2 (x))
=p.
Remark F2. We define kp -interval as the p-CI for constructed as (I (z), us
(x)) for each x E Q
where
L 2 (x) = min
(1 - 20 + 221,
||s- Bll
(23)
U (x)= = arg max
(1 - 20 + 221-
The value 13 is the squared root of the upper (1 - p)/2 critical value of the chi-squared distribution with 2 degrees of freedom.
Proef: According to Theorem
Pr(fL 22 (z) < u(() < 23 Yo v
Pr([110' Z112 kx²)))
(24)
=p.
The left side of Figure F.1 provides an indicative example of a kg -interval and kz -interval, with = 95%
and 8" = (0, 0). The results are obtained using only one sample Z = (Zo, 71), which can be visualized by considering that (0) Zo and u1()) Z1. Further, the empirical confidence level of an interval (T, U) can be estimated as == $ > 1( LTz (ke) < u(ke) < Uz, (ke)}) with > 0. That is, we calculate proportion of times that such a construction of the interval contains the true values on [1/e]e]- According to law large numbers, the empirical confidence level converges with increasing J. probability, to the true confidence level p. The right side of Figure F.T shows that the empirical confidence level curve of the b33" -interval (or the k. -interval) converges a value above (or bclow) 95%. This result implies that the k.x2 -interval is valid 95%-CI -I for while the k.Y -interval is not.
16

---

Veriaa
Doed
0375
0200
66 2.00 208
1899
ulo
Figure F.1: Left: An indicative example of ky inteerall and a kx iineerval, with p 95% and 0 = 0). The results are obtained using only one sample Z = (Zo. Z1). which can be visualized by considering that us (0) Zu and (1) Z1. Right: Empirical confidence level on versus the number of samples J.
G Proofs
G.1 Proof Theorem 4.1
If E ž,, then L(x) = min.ež, u(z; 2) I ua(x;)) & max,c2, u(x;2) = Uxx). for all x Ω. Further, according to the definition of the proper solver, we have u1(2; Σ) - S u(x) < ua(2; 2) + for all E Q, and thus Lxx) s u(x) < (2x) + m, for all x Ω. Finally, can conclude that Pr(L(x) 7 S u(x) < (xx) + Vo E tt 2 Pr(z E 2,)2
H Algorithms
The developed methods summarized in Algorithms and
Algorithm PICProp: Physics-Infooe Informed Confidence Propagation (Gradient-bare based)
Input: (p-CC for 2), query point z, inner-steps schedule [N& Ke steps Ko.
k= learning rates a, B, outer
Initialize 2g E žp. 0g.
for A 1,..., Ko do
for n = 1,... N₁ do
8, 0 V Cpinninnnnn. (fn- 2k-1)
end for
i 0g ON
if lower bound sought then
2k 2k-1 BV M(== 1) (I)
else
=k-1 + BV 3 (x:
end if
end for return (2(2x.)
I Bi-level Optimization
A bi-level optimization problem can be formulated as
min J(0 (0). +) s.t. (
(())
min [(0, x),
(25)
17

---

Algorithm 2 EffiPIPPrpp: Efficient Physics-Informed Confidence Propagation
Input: p(agg, X. Randomly sample sr K ~ p(xg))
Obtain (x)(()))) ((= (xf)) with Algorithm for 1, ... K.
Initialize sh. while not done do
/ db v/ Aa ;p) according to Equation 7
end while return Up
where 0, á are the parameters the hyperparameters, whereas and are lower-level objective and upper-level objective. Many problems like meta-learning, hyperparameter optimization, data poisoning, and auxiliary learning can be naturally formulated in bi-level optimization way.
The hardest part of bi-level optimization is that 0" (() is not availahle in closed-form in many scenarios, particularly with empirical lower-level objective that requires multiple steps of gradient descent. In such cases, the explicit gradient of high-level objective with respect hyperparameter intractable.
According the chain rule.
(26)
The most challenging part is the approximation of the implicit gradient v000 (o), for which many advanced methods have been proposed studied recent years [15). this paper. three methods referred as Reverse. approximate implicit differentiation with Neumann Series (AID-NS), and approximate implicit differentiaion with conjugate gradient (AID-CG) are evaluated on the confidence interval identification tasks.
1.1 Reverse
Reverse is an iterative differentiation bbasiatto----dddddddee method to approximate the hypergradient by unrolling the inner optimization 01:K trajectory in reverse way. The implicit gradient can be approximated by
K
00:+1
00
Vo0' Vo)xx()))
(27)
00:
1=++1
where 0++1 = 0k - aVol(Os <) and a is pre-detined inner-learning rate.
1.2 Approximate Implicit Differentiation (AID)
Another family of methods to approximate the hypergradient is Implicit Differentiation. where the implicit function theorem (IFT) used evaluate va0'(d If second-order derivable and strong convex, we have
Va0" (0) = - -
(28)
Then approximate the gradients of using
=)
(29)
The most challenging part is the approximation of the inverse Hessian vector product (inv-hyp) A 1 where Raxd is a given symmetric positive matrix and p R is given vector.
AID-NS 19 is method based on the Neumann Series to approximate the inv-hvp as
K
A A p
-
Σ XA)* pm
(30)
k=0
R=0
where À is hyperparameter that ensures convergence and K is hyperparameter for precision control.
AID-CG [12] approximates the hypergradient by solving the system of homogeneous linear equations Az with conjugate gradient.
18

---

J Details regarding implementation, hyperparamelers, query points, and computational times
Implementation and hyperparameters J.1
For each PDE, we use multilayer perceptron (MLP) for the PINN model and another MLP of the same architecture for EfoPICPrrp. Query points collected in grid manner. For each query point, PICProp starts with standard PINN training with a randomly initialized virtual boundary condition. We refer this as warmup phase. Subsequently. PICProp alternates between performing a number inner steps to update the PINN model and computing the hyper-gradien for meta step to update the virtual boundary condition. The inner/outer learning rate selected from {0.001, 0.01, 0.1]. The inner/outer optimizer selected from SGD and Adam. The hyper-gradient method is selected from Neumann Series (NS), Conjugate Gradient (CG), Reverse; warmup/inner/mett steps are tuned balance stahility
sec more details Appendix | The numbers of bi-level training and the computational cost. For the pedagogical example, as six points are too few for train- validation splitting, we manually set À 1 for EfiPiCPrrp. In the rest of the examples, split 10% of the training data the validation set select the best À from {0.0, 0.25, 0.5, 0.75, 1.0). The hyperparameters are summarized Table I.
Table J.1: Summary of implementation details, and utilized hyperparameters and architectures in the experiments.
Pedagogical
Poisson
Burgers
MLP arch.
32 x
20 8
20 x
# BC points
2
40
256
# query points
6
11 11
41
inner optim.
Adam
Adam
Adam
inner
0.001
0.001
0.001
meta optim.
SGD
SGD
SGD
meta Ir
0.01
0.001
0.001
hypergrad
NS
Reverse
Reverse
warmup steps
2,000
4,000
20,000
# inner steps
500
20
100
# meta steps
50
200
500
X
1.0
1.0
0.0
J.2 Query Points
The query points are summarized in Figure J.1
(a) Poisson 2D: Query Points
(b) Burgers: Query Points
1.0
1.0
5
0.75
0.50
X
X
4
0.0
3
0.0
0.00
0.25
-055
2
-0.5
0.50
-0.55
X
X
1
-1.0
-1.0
-1.0
Figure J.1: Visualization of the query points on top of the exact solutions of Equations II and 13 (a) 2-D Poisson equation: 121 equally spaced points grid. (b) Burgers equation: 246 equally spaced points on a grid.
19

---

J.3
Computational times
our computational experiments, methods were run serially and thus, the runtimes shown below are much larger than the required ones for using these methods in practice. For this reason, we have also included "pseudo-paralle" runtime that shows the runtimes if parallel processing had been used instead. The runtimes for the examples of Section are summarized in Table
Table J.2: Computational time for different experiments.
Experiments
ES
PICProp
SimPICProp
EffiPICProp
Pedagogical
42 hrs for 5,000
1.7 hrs for
16 mins (40s
18 mins (40s
example
trials
queries
PICProp queries)
lambda + 6 PICProp queries)
(CPU)
2D equation (GPU)
Poisson
14 hrs for 1,000 trials
400s per query
14 hrs (180s per lambda + 121 PICProp queries), 10 mins for pseudo-parallel
70 hrs for 5,000 trials
ID
time-
347 hrs for 2,500 trials
5,000s per query
14 hrs (2000s per lambda 246 PICProp queries), 2 hrs for pseudo- parallel
dependent
Burgers equa- tion (GPU)
694 hrs for 5,000 trials
K Additional Experimental Results
K.1
Pedagogical example
K.1.1
Learning Curves
provides meta-leaming curves corresponding predictions Figure 2 (main text). provides the results all considered meta-learning methods and ES on six query points. We conclude that the combination of Neumann Series (NS) and SGD stably converges to the reference result of ES.
Figure K.
Figure 2
K.1.2
Figure K.3
provides comparison among PICProp and PrfiPCCPopp with various As. Notice that EffiPCPrrpp
with X is referred SimPICProp in Figure 3 The mean squared errors (MSF) hetween PlCProp CI and
EMPICCPrP CI summarized in Appendix K.1.4 We conclude that EffiPIPrrpp with 1.0 is the best approximation PICProp Cl.
K.1.3 SimPICProp EMiPICProp
Figure K.4 and Figure K.5 provide comparison between SimPICProp EMPICCroP T and Hoeffding boundary CI, respectively.
K.1.4 Validity of EMAPICPPrp Cls
Although has been already presented in Figure 2 that PICProp and EffiPICProp predictions well match the ES results, we conduct the following experiment to further verify the validity of EfriIICPrpp Cl predictions.
We independently sampled 500 boundary conditions and propagate the corresponding x² CIs using EMPICCPrp to construct CIs for the entire domain. According to the definition of joint CIs (see Definition 3.1 instance of CI is valid if it contains the solution u(r)) for any E [-1. 1). In this regard, we uniformly sampled 101 data points from [-1, 1] and regard an instance of as valid the solution was enclosed within the interval for these data points. Further, we run vanilla PINNs with clean boundary data and use the max absolute error among these 101 data points as As result, 486 out of these 500 197.2% > 95%) instances of CI are valid, which empirically indicates that EttiPICProp CI predicitons valid and conservative.
20

---

2 -02
-0.3
340
945-
Opicessinn
0.0
16:0
-425
-1.00
-0.6
-020
-1.2
-0.78
10
20
20
40
50
16
20
30
46
10
23
30
mada
stgs
0.650
0.10
0.625
1.20
0.45
0.800
8.75
0.00
20
PiChrpp (NS S0/01
Figure K.1: Pedagogical example: The learning curve of PICProp with Neumann Series (NS) and SGD six query points. Comparisons with exhaustive search results for solving the bi-level optimization problem for each query point are also provided.
-10
-2
-0.5
x- 0.2
2-0.8
2= 1.8
1.05
0.32
0.825
1.50
0.05
0.600
0.676
0.65
0.00
Figure K.2: Pedagogical example: Various combinations of grreietgradeett methods with meta- optimizers and comparisons with exhaustive search results for solving the bi-level optimization problem for each query point.
À MSE
0.0
0.25
0.5
0.75
1.0
1.4E3
2.2E-4
2.2E.4
1.7E4
18E-5
Table K.1: Mean Squared Errors between PICProp CI and EffPICPrrp CL. EffiPIPProp CI with 1.0 the best approximation of PICProp CI.
K.2 2-D Poisson equation
Figure/K. K.6lincludes slices along the axis and accompanies Figure (main text), which includes slices along the axis. Heatmaps with Cls over the entire domain presented Figure K.7
21

---

Chi-squared
Chi-squared (zoomed
1.0
RDY LOROO -ut
0.5
0.0
hchhp
DTFiCep 4-
1-12
-1.0
Figure K.3: Pedagogical example, x2: Comparison among PICProp and EffiPICProp with various
Às. Notice that EftiPICProp with À = 0 is referred to as SimPICProp in Figure
T-squared
T-squared (zoomed in)
Friml ghfichop
1.2
1.0
18728
1.0
0.8
0.6
CETIChOP
0.4
-1.0
-0.5
0.0
0.5
1.0
0.2
0.3
0.4
0.5
0.6
Figure K.4: Pedagogical example, 72: Although the results of both SimPICProp and EfriPIPProp
match PICProp results on the query points, EffiPICProp generalizes better on unseen locations.
K.3 Burgers equation
Figures
K.8 & K.9
provide the meta-learning curves corresponding to the CI predictions of Figure 3 (main text).
L Comparison with Bayesian PINNs
In this section, we compare PICProp with Bayesian PINNs (BPINNs) 3 for the problem of Section 5.1 For BPINNs, the inference method for posterior estimation is Hamiltonian Monte Carlo (HMC) 4. Although provided here completeness, we believe that such comparison can also misleading because of the many differences between the two methods. The first difference pertains to the fact that PICProp does not consider parameter uncertainty. Further, the assumptions for these two methods are different: BPINN requires heavy assumptions about the data likelihood, for instance, while PICProp only requires the noise to be mean-zero. Furthermore, PICProp produces joint CIs for the whole domain, while the BPINN intervals marginal. Finally, although supported by sound Bayesian justification, BPINNs are not supported by a strict confidence guarantee for the produced intervals, in contrast to Theorem 4.1 that supports PICProp.
We summarize and contrast the two approaches in the following: BPINNs assume a NN parameter prior, and subsequently approximate and sample the parameter posterior to obtain samples from the solution posterior. This is done in order to fit solution posterior distribution and subsequently compute its statistical intervals. In contrast, the problem of UQ is rethought ab initio our paper. PICProp first considers either CIs given at the boundarics or data collected the boundaries and estimates the corresponding Cls. Subsequently, the boundary CIs are propagated to the interior domain with mild assumptions and with the theoretical guarantee that the obtained CIs are indeed CIs according to the definition.
From the results summarized in L.1 observe from Figure L. that the PICProp CIs more conservative than the BPINN. From design perspective, conservative uncertainty estimation is the preferred route for critical risk fisksvensitive applications 7 although can be more costly. In contrast, most widely used UQ methods. such as deep ensembles [88] and Monte-Care dropout (28 10, are known to be overconfident in practice 7 Finally, note that the BPINN results of Figure L.] were obtained by using the software NeuralUQ developed in 88
22

---

Hoeffding
Hoeffding (zoomed in)
SEOP Chop
1.0
PR LEBSCh03
0.0
-0.5
0.4
sarch:
-1.0
-0.5
0.0
0.2
Figure K.5: Pedagogical example, Hoeffding: Although the results of both SimPICProp and EffiPICProp match PICProp results on the query points, EftiPICProp generalizes better on unseen locations.
[a.1] 0.5
[a.2] y Izoomed in)
y
0021 y 0.0 (zoomed
35
2.5
23
20-
22-
71-
15
28-
1.8
[c.1] 0.5
(c.2) 0.5 (2oamed in)
(d.1) y- -1.0
(d.2) 1.0 (zoomed in)
0.5
4.-
4.0
-25
Exact
EffiCChopp (121 points)
(1000 points)
15000 points)
Figure K.6: 2-D Poisson equation: Slices of CI predictions for y E 0.5, 0.5, 1.0).
20
Figure K.7: 2-D Poisson equation: Heatmaps of upper and lower CIs predicted by the exhaustive search (ES) and EffiPICProp.
50
P0003- 200 1380
f 0.0
(b) f 0.25
(c) f = 0.5
(d) f = 0.75
Figure K.8: Burgers equation: Meta-learning curves for 0 demonstrating the stability of our method.
23

---

(a) (
0.0
(b) ( 0.25
(c) [ 0.5
(d) ( 0.75
Figure K.9: Burgers equation: Meta-learning curves for = 0.2 demonstrating the stability of our method.
Chi-squared
sintmx
1.0
EMRPICProp
BPINNS+HMC
sample
00.5
-1.0
-0.25
0.00
0.75
Figure L.1: Pedagogical example: The comparison between intervals produced by EMPICProP and BPINNs [34|+MC [4]. Noise on the boundary is Gaussian with known standard deviation 0.05.
24


Reviews:
[{"review_id": "Sfh2Jy8vAX", "reviewer": "Reviewer_H2bG", "summary": "The paper introduces a PICProp, a method to propagate uncertainty through spatiotemporal nonlinear PDE systems. PICProp (1) is more valid and makes less assumptions compared to NN based UQ methods such as dropout and deep ensembles, and (2) can be modified to be more computationally efficient by learning a meta-model. The paper provides extensive computational to show the effectiveness of PICProp.", "questions": "See weaknesses.", "limitations": "N/A", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The methodological section is very clear and easy to follow. PICProp is a principled, concise, and elegant method for the proposed problem setup. \n- I really enjoy the writing and organization in the experiment section as well. The pedagogical example is illustrative of how the algorithm works for different types of CIs, and the following sections showcase the effectiveness of the method in 2d and in PDEs with steep gradients. ", "weaknesses": "1. The aper mentions that PICProp produces conservative CIs that are \"clearly preferred in critical and risk sensitive applications\". I don't think this statement is necessarily true. In many cases, one wants the CIs to be valid and **informative**. You can have very large CIs that are trivially valid but too wide to have any meaning in decision-making. (In general, wide CI = safety is neither a correct or scientifically sound statement, don't use it in your papers. Line 285 is an example.) There exist metrics such as Winkler Score and interval scores [1] that principally evaluate the quality of a CI; I recommend the authors to use them. The authors did mention this shortcoming in the conclusion section and say that they will study it in a future work - My opinion is that these studies need to be a part of *this* work for the method to be convincing and usable for the community.\n\n[1] Gneiting, Tilmann, and Adrian E. Raftery. \"Strictly proper scoring rules, prediction, and estimation.\" Journal of the American statistical Association 102.477 (2007): 359-378.\n\n2. The marginal probability framing of CIs (Theorem 4.1) is very close to the guarantees provided by frequentist UQ methods (jackknife and conformal prediction, for example). These are also \"standard\", distribution-free, and widely used UQ methods; it's important to mention and compare to them. Comparing to them will likely help you prove your case of why it's important to have CIs be physically informed. \n\n3. This might be more from of a machine learning community perspective, but I'm curious to see PICProp's performance on a real dataset instead of just on synthetic experiments. The applications authors mentioned in appendix A are all good candidates.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful feedback. Here, we provide answers to your concerns point by point.\n\n**To Weakness 1:** We would like to clarify that none of the statements in our paper are intended to give readers the impression that \"wider CI = safety.\" In fact, in practical applications, constructing confidence intervals involves a trade-off between practicability and safety. \nTypically, we aim for informative and valid confidence intervals to aid in making more accurate decisions. However, in **risk-sensitive** applications, considering the significant harm of non-valid confidence intervals, constructing conservative confidence intervals by sacrificing some information to ensure the validity might be preferable.\n\nGiven that commonly used uncertainty quantification (UQ) methods, including deep ensembles and MC dropout, tend to be over-confident (too narrow to be valid) in practical applications, we emphasize the conservatism of PICProp as a distinguishing property from traditional approaches. Users can select accordingly based on the specific requirements of their practical applications.\n\nRegarding revision, we will follow the reviewer's suggestions to (1) enhance the comprehensiveness of the relevant discussions in the introduction and (2) revise other parts of the manuscript where there might be statements that inadvertently suggest \"Wide CI = safety,\" including but not limited to line 285, line 309 and the caption of Figure 5.\n\nWe agree with the reviewer that universal metrics for evaluating confidence intervals will help to enhance the comprehensiveness of this work. However, interval scores defined in [1] are for **prediction intervals** rather than confidence intervals. To the best of the authors' knowledge, such a metric for confidence intervals has not been well studied yet.\n\nWe are open to incorporating general metrics to evaluate confidence intervals in future research. For example, a straight-forward generalization of interval scores to confidence intervals may be $S(L, U, z) = U(z) - L(z) + \\frac{\\alpha}{2}(L(z) - u(x))\\mathbb{I}\\{u(x) < L(z)\\} + \\frac{\\alpha}{2}(u(x) - U(z))\\mathbb{I}\\{u(x) > U(z)\\}$, where $z$ denotes an instance of boundary data, and $(L(x), U(x))$ denotes the constructed CI for the solution $u(x)$ at $x$. But still, these potential metrics need to be thoroughly studied before being applied to our work.\n\n**To Weakness 2:** We would like to highlight that the problem studied in this work is to construct **confidence intervals** for the solution of PDEs, which is constant as we are considering deterministic systems. Jackknife and conformal prediction mentioned by the reviewer, are approaches to construct \\textbf{prediction intervals} for the next sample, which is indeed a random variable. Thus, although considered as standard UQ methods, they cannot be directly applied to construct CIs.\n\nWe also would like to direct the reviewer to Appendix B, C, and L for more detailed discussions about the difference between PICProp and traditional UQ methods (e.g. Bayesian method, MC simulation), especially the differences in the problems they address and the properties of the probabilistic intervals they construct. Overall, the results obtained by these methods can be significantly different both in terms of interpretation and quantitatively. \nAs a result, we decided to omit such a comparison for avoiding any potential confusion of the readers. \n\n**To Weakness 3:** The dimensionality of PDE/ODE problems is typically N+1, where N is the number of spatial dimensions and 1 accounts for the temporal dimension. In this regard, the problems that we considered in this paper involve standard representative equations met in many real-world problems including heat transfer, fluid dynamics, flows through porous media, and even epidemiology [1-2] and N is in many practical cases equal to 1 or 2. The problems that we selected to present in this paper are similar to or the same as the standard problems solved by the computational science community for developing novel solution methods (both deterministic [3] and considering uncertainty [4]).\n\n[1] Meng, X., Babaee, H. and Karniadakis, G.E., 2021. Multi-fidelity Bayesian neural networks: Algorithms and applications. Journal of Computational Physics, 438, p.110361.\n\n[2] Zou, Z., Meng, X., Psaros, A.F. and Karniadakis, G.E., 2022. NeuralUQ: A comprehensive library for uncertainty quantification in neural differential equations and operators. arXiv preprint arXiv:2208.11866.\n\n[3] Psaros, A.F., Kawaguchi, K. and Karniadakis, G.E., 2022. Meta-learning PINN loss functions. Journal of computational physics, 458, p.111121.\n\n[4] Yang, Y. and Perdikaris, P., 2019. Adversarial uncertainty quantification in physics-informed neural networks. Journal of Computational Physics, 394, pp.136-152.\n"}, {"review_id": "ZBVFY8Uspa", "reviewer": "Reviewer_Lv8K", "summary": "The paper discusses the limitations of current approaches for uncertainty quantification in deep learning and physics-informed learning, such as the need for strong assumptions regarding data likelihood and the computational cost associated with approximating the posterior. The paper introduces a novel problem of confidence interval estimation for deterministic partial differential equations and proposes a method called Physics-Informed Confidence Propagation (PICProp) based on bi-level optimization to compute a valid confidence interval without making heavy assumptions. The paper provides a theorem regarding the validity of the method and computational experiments focusing on physics-informed learning.", "questions": "Similar to the weakness, I have several concerns about the claims and experiments demonstration \n\n1. The experiments can be solved by some baseline methods as the author mentioned because the cost is acceptable. If the authors would like to show the real power of this method, they may consider some real-world problems, either from PDE or simulations. \n\n2. The computational cost is the major barrier to the UQ problem. The method leverages bi-level optimization, how about the cost compared with the traditional methods, e.g., MC-based approaches, specifically in complex cases? For example, high-fidelity simulations, CDF or FEM, can you handle the optimization efficiently without any approximations?  How about the computational cost? \n\n3. I only saw the analytical experiments and have been thinking about how to apply this method to solve real-world applications, is that useful or helpful to handle complex cases in AI for science problems? \n\n----\nAfter rebuttal\n\nI posed several follow-up questions but unfortunately, I did not receive the authors' positive responses although they made efforts to argue. As a result, I incline to keep my original score.", "limitations": "Did not see a clear discussion about the limitations.", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The paper is simple and clear to follow\n2. The paper is well structured \n3. The paper aims at solving an important problem involving UQ in PINN\n4. The theoretical analysis is helpful to better understand the benefits of the method", "weaknesses": "1. The paper claims efficiency but did not demonstrate the cost saving with enough supports\n2. The scope is narrowed to PDE but did not see practical applications, e.g., large-scale cases \n3. The experiments only focus on toy examples that did not show the strength of this work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful feedback. Here, we provide answers to your concerns point by point.\n\n**To Q1:** We assume that the reviewers are referring to the exhaustive search (ES) method as the baseline, of which the cost is acceptable. \n\nWe would like to clarify that the cost of ES is acceptable only for the pedagogical example (2 boundary data points). \nWhile acceptable, the cost of ES (42 hrs for 5k trials) is several times higher than that of EffiPICProp (less than 20 mins) to achieve a similar approximation of the solution of Eq(4).\nAs the number of data points increases, e.g. 2D Poisson equation (40 boundary data points) and Burgers equation (256 boundary points), the cost is not acceptable anymore. \nEven with a cost budget several times larger than that of PICProp, ES still struggles to well approximate the solution of Eq(4) accurately. One possible reason is that as the number of data points increases, the search space grows exponentially, leading to a corresponding exponential increase in search costs.\nIn this regard, the authors conclude that PICProp is more powerful and scalable compared to ES, in terms of solving Eq(4).\n\nIf the reviewer meant traditional UQ approaches as baselines, we would like to direct the reviewer to Appendix B, C, and L for more detailed discussions about the difference between PICProp and traditional UQ methods (e.g. Bayesian method, MC simulation), especially the differences in the problems they address and the properties of the probabilistic intervals they construct. Overall, the results obtained by these methods can be significantly different both in terms of interpretation and quantitatively. In this regard, these approaches cannot serve as baselines for comparison with our approach.\n\nMoreover, the problems that we considered in this paper actually met in many real-world problems. Please also refer to the answer to Q3.\n\n**To Q2:** We itemize our response to this question as follows: \n\n(a) PINNs as a deep learning approach for solving PDEs are an alternative to conventional CDF or FEM simulations.\nAs the focus of this paper is to use PINNs as the underlying numerical solver of the PDE to construct CIs, the comparison between PINNs with traditional simulation approaches is outside the scope of this paper and has been discussed in [1].\n\n[1] Karniadakis G E, Kevrekidis I G, Lu L, et al. Physics-informed machine learning[J]. Nature Reviews Physics, 2021, 3(6): 422-440.\n\n(b) A seemingly feasible alternative to our approach is Monte Carlo simulation using multiple FEM simulations. However, it has the same limitations as Monte Carlo simulation using PINNs as the PDE solver, which are discussed in Appendix C.\n\n(c) We would like to emphasize that there are many differences between PICProp and traditional UQ methods (e.g. Bayesian method, MC simulation), as discussed in Appendix B, C, and L.\nOverall, the results obtained by these methods can be significantly different both in terms of interpretation and quantitatively. \nAs a result, we decided to omit such a comparison for avoiding any potential confusion of the readers. \n\n(d) Gradient-based approaches used in this work for bi-level optimization, are all approximate methods. \nPlease refer to Appendix I for more details on approximating the implicit gradient.\nIn terms of computational cost, all three methods for approximating implicit gradient used in this paper provide hyperparameters to balance computational cost and accuracy.\nThe performance of these methods depends on the specific problem and is difficult to analyze.\nFor the problems studied in this paper, we summarize the empirical computational cost in Table J.2.\nCompared to the search-based method, bi-level optimization-based methods\nare more efficient.\n\n**To Q3:** The dimensionality of PDE/ODE problems is typically N+1, where N is the number of spatial dimensions and 1 accounts for the temporal dimension. In this regard, the problems that we considered in this paper involve standard representative equations met in many real-world problems including heat transfer, fluid dynamics, flows through porous media, and even epidemiology [1-2] and N is in many practical cases equal to 1 or 2. The problems that we selected to present in this paper are similar to or the same as the standard problems solved by the computational science community for developing novel solution methods (both deterministic [3] and considering uncertainty [4]).\n\n[1] Meng, X., Babaee, H. and Karniadakis, G.E., 2021. Multi-fidelity Bayesian neural networks: Algorithms and applications. Journal of Computational Physics, 438, p.110361.\n\n[2] Zou, Z., Meng, X., Psaros, A.F. and Karniadakis, G.E., 2022. NeuralUQ: A comprehensive library for uncertainty quantification in neural differential equations and operators. arXiv preprint arXiv:2208.11866.\n\n[3] Psaros, A.F., Kawaguchi, K. and Karniadakis, G.E., 2022. Meta-learning PINN loss functions. Journal of computational physics, 458, p.111121.\n\n[4] Yang, Y. and Perdikaris, P., 2019. Adversarial uncertainty quantification in physics-informed neural networks. Journal of Computational Physics, 394, pp.136-152."}, {"review_id": "fAFpJM2GDb", "reviewer": "Reviewer_Jho9", "summary": "Authors present a method to propagate confidence intervals from some locations to the entire domain. In specific, they consider the problem of propagating uncertainty / confidence internals at the boundary, in partial differential equations based systems, to the rest of the domain. The proposed method PICProp is a bi-level optimization method. Authors present theoretical guarantees for PICProp and then develop a computationally efficient version of the PICProp under some simplifications. The authors show the efficacy of the method via some computational experiments.", "questions": "It would be great if authors can comment about points mentioned in \"Weakness\" section", "limitations": "Please refer to \"Weakness\" section. \n\nEven though the points mentioned in the weakness section can be thought of as limitations, I don't find them to be major.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper tackles an important problem of estimating and propagating uncertainty in physics-informed ML systems. The proposed methods are the first methods (to the best of my knowledge) to develop confidence interval propagation method for PDE systems.\n\nThe authors present a method, Physics-Informed Confidence Propagation (PICProp), to propagate confidence intervals from any boundary points in the domain to any other points in the domain. They show theoretical guarantees which demonstrate that PICProp propagates confidence intervals with high probability.  The theoretical analysis offers insights into the inner workings of PICProp.\n\nThen authors identify that PICProp is computationally expensive and might not be applicable beyond small domains. To address this they develop an efficient version of PICProp known as EffiPICProp which makes confidence interval propagation much more computationally viable. It is great that authors through about the practicality of the algorithm and developed this. \n\nAuthors test their method on computational experiments. In spite of being small scale experiments, they offer intuition which is important for transferring these methods to other domains / problems. It is great that the authors have actually documented the computational time for running different methods on these problems.\n\n", "weaknesses": "All the computational experiments seems to be conducted on 2D domains. It is not clear how the computational costs scale once we go to higher dimensions and if EffiPICProp can still be computationally feasible. It might be useful to discuss about this.\n\nThis work deals with deterministic systems and often most real-life applications involve stochastic systems (for example, due to inaccurate modeling of the systems). It might be useful to discuss about how the proposed methods can or cannot transfer to stochastic systems.\n\nAs of now, the proposed method only propagates uncertainty from boundary points. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful and overall positive feedback. Here, we provide answers to your concerns point by point.\n\n**To Weakness 1:** The dimensionality of PDE/ODE problems is typically N+1, where N is the number of spatial dimensions and 1 accounts for the temporal dimension. In this regard, the problems that we considered in this paper involve standard representative equations met in many real-world problems including heat transfer, fluid dynamics, flows through porous media, and even epidemiology [1-2] and N is in many practical cases equal to 1 or 2. The problems that we selected to present in this paper are similar to or the same as the standard problems solved by the computational science community for developing novel solution methods (both deterministic [3] and considering uncertainty [4]).\n\n[1] Meng, X., Babaee, H. and Karniadakis, G.E., 2021. Multi-fidelity Bayesian neural networks: Algorithms and applications. Journal of Computational Physics, 438, p.110361.\n\n[2] Zou, Z., Meng, X., Psaros, A.F. and Karniadakis, G.E., 2022. NeuralUQ: A comprehensive library for uncertainty quantification in neural differential equations and operators. arXiv preprint arXiv:2208.11866.\n\n[3] Psaros, A.F., Kawaguchi, K. and Karniadakis, G.E., 2022. Meta-learning PINN loss functions. Journal of computational physics, 458, p.111121.\n\n[4] Yang, Y. and Perdikaris, P., 2019. Adversarial uncertainty quantification in physics-informed neural networks. Journal of Computational Physics, 394, pp.136-152.\n\n**To Weakness 2:** We thank the reviewer for the constructive suggestion.\nThe solution of a stochastic system is itself stochastic.\nIn this regard, confidence intervals are no longer applicable, as confidence intervals are defined for unknown constant statistics rather than random variables.\nAlthough the problem formulation of uncertainty quantification for stochastic systems is totally different, the idea to propagate the uncertainty might be adapted.\nFor example, if we can construct CIs for the force terms, we can propagate the uncertainty from force terms to the solutions by including the force terms into the bi-level optimization in the exact same way as the boundary conditions. \nWe leave the more detailed problem formulation, theories, and algorithms for future work.\n\n**To Weakness 3:** PICProp can be easily generalized to propagate uncertainty from inner data points or collocation data points, once CIs/data to construct CIs at these points are available. We plan to work on it in the future."}, {"review_id": "JFZallwSPY", "reviewer": "Reviewer_JUgH", "summary": "The paper focuses on a novel algorithm for confidence interval estimation for deterministic partial differential equations. It provides leverages bi-level optimization to compute a valid CI, that is supported by both theoretical and numerical investigations, focusing on physics-informed learning.", "questions": "1. Please complete table J.2 and the show the other comparisons in the figures 3,4,5 of SimPICProp and PICProp to ES and EffiPICProp.\n2. In appendix C, you should mentioned the well-known result of Jensen's inequality as the figures are an illustration of this.\n3. Please consider discussing the Bayesian PINNs in the main text.\n", "limitations": "Yes", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is very well written, with clearly structured sections and appendices. It provides a lot of technical detail and manages to effectively balance what to be included in the main text and what to be pushed to the appendices in order to not distract. I particularly enjoyed reading this work as in blends theoretical aspects with computational ones, and has direct and immediately applicability to the computational science world, particularly error estimations.\n\nThe authors have performed a number of benchmarks (of increasing complexity) for the methodology, and have clearly showed its advantages as well as observed (\"We conclude that EffiPICProp, although computationally cheaper, in most cases provides wider [...] CIs, as compared to ES with varying numbers of trials\") or potential drawbacks (\"Although not encountered in the experiments, a potential limitation of our method is that the obtained CIs can be loose in some cases, even when propagating a tight CI from the boundaries, due to the max/min operations in the bi-level optimization of Equation 3\").\n\nOverall, I believed the presentation of the work is of high quality and has impact implications to its directed target.", "weaknesses": "A weakness of the paper is the fact that some benchmark results are incomplete. Notably SimPICProp in the pedagogical example is only shown for Chi-squared, not for T-squared or Hoeffding. Also SimPICProp and PICProp are not shown for the 2D Poisson equation or the 1D time-dependent Burgers equation.\n\nLastly, while as the authors say \"The aim of this paper is to propose a principled and theoretically guaranteed solution approach and thus, we refrain from applying such ad hoc solutions in the computational experiments and plan to consider parameter uncertainty rigorously in a future work. Nevertheless, we provide in Appendix L a comparison with Bayesian PINNs for completeness.\", the figure in Appendix L suggest a potential better performance of the Bayesian PINNs compared to EffiPICProp. For objectivity reasons, I would suggest moving some of that discussion into the main text.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful and overall positive feedback. Here, we provide answers to your concerns point by point.\n\n**To Weakness 1 and Q1**: We thank the reviewer for pointing out the absence of some results. We itemize our response to this question as follows: \n\n(a) We will show SimPICProp results in the pedagogical example for T-squared or Hoeffding CIs in the appendix of the revised manuscript.\n\n(b) SimPICProp is a special case of EffiPICProp when $\\lambda$ equals 0. \nWe use the term \"SimPICProp\" for the naive regression approach to distinguish it from the meta-learning approach, EffiPICProp ($\\lambda > 0$). \nFor 2D Poisson equation and 1D time-dependent Burgers equation, we split $10\\%$ of the training data as the validation set to select the best $\\lambda$ from $\\{0.0, 0.25, 0.5, 0.75, 1.0\\}$ (see line 606).\nNotice that $\\lambda=0$, i.e. SimPICProp is in the candidate set.\nIn this regard, SimPICProp results and costs are not separately reported, as SimPICProp is a special case of the principled EffiPICProp.\nFor revision, we will include EffiPICProp results with different $\\lambda$s, because we believe such results will help the reader to better understand the power of meta-learning.\n\n(c) The reason we omit the result of (brute-force) PICProp for 2D PDEs is that the cost is too expensive. For example, the grid we use for plotting the result for Burgers equation is $256 \\times 100$, which requires $25600$ queries. Each query takes 5000s, as reported in Table J.2. The total computational cost would be 4 GPU years, while EffiPICProp only takes 14 GPU hours.\n\n**To Weakness 2 and Q3:** We thank the reviewer for the revision suggestion. We will accordingly move some of that discussion into the main text in the revised manuscript.\n\n**To Q2:** We thank the reviewer for the constructive suggestion. We will mention Jensen's inequality in the revised manuscript."}, {"review_id": "sK2u0v4yhz", "reviewer": "Reviewer_5z2U", "summary": "The paper presents methods to estimate confidence sequences for the solution of a differential equation given noisy boundary conditions. Various methods are proposed of different complexity and the confidence sequences are illustrated for the solution of three different problems. ", "questions": "In eqn. (5) it seems that the lower and upper boundary data points should be derived by finding the points that  minimized u_{\\hat{\\theta}}(x) at x = x_q^k rather than for all x.  Is that correct?\n\nFour lines above eqn. (4), should z minimize (rather than maximize) the corresponding value to determine L(x)?\n\nThird to last last line of section 4, SimPICProp is not defined in the main paper.  It may be better to state that SIMPICProp is the version of EffPICProp where \\lambda = 0. \n\nThe paper use the term confidence interval for Def. 3.1.  Based on the work of Ramdas, I believe this is a confidence sequence. The authors should verify if this is indeed the case.", "limitations": "No negative societal impact as this is foundational work that is not considering a specific application or using an datasets.  The actual differential equations being solved are toy problems.\n\nLimitations of the method were addressed by the review in terms of capturing the accuracy of the neural solution.  This limitation was not discussed in the paper. ", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The overall problem to determine the confidence sequence is laid out.  Under the given assumptions, the problem is well formed and an approximation of the exact solution via exhaustive search is clear.  The faster methods that incorporate meta-learning techniques and surrogate solutions for the lower and upper sequences are innovative and clever. ", "weaknesses": "The main assumption of the paper is that the accuracy of the solution provided by the machine learning model is either perfect (\\eta = 0) or \\eta is known and fixed when the boundary conditions are given perfectly.  Even with clean boundary data, the accuracy of the model will depend of the number of internal and boundary training points, and these are not incorporated in the computation of the confidence sequences.\n\nWhile it might be OK to consider the exhaustive search as the gold standard, the examples do not confirm whether or not the derived confidence sequences are indeed satisfying Definition 3.1. In light of the first comment, this should be validated. \n\nAlgorithm 1 and 2 are part of the main contribution of the paper.  However, they are tucked away in the supplementary appendix. As such, a reader if forced to go to the supplementary material to even understand the gist of the paper. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your thoughtful feedback. We will address your comments point by point below, and we hope that our answers can ease your concern.\n\n**To Weakness 1:** In terms of not considering epistemic uncertainty, the reviewer is right and we have already clarified in the paper that we consider only proper solvers (see Definition 3.3) and thus, the focus is placed on the propagation of aleatoric uncertainty. \nMore specifically, epistemic uncertainty is simplified to be known and fixed as $\\eta$.\nAs for the specific quantification method for $\\eta$, it is beyond the scope of this article's discussion. \nMoreover, any quantification method for $\\eta$ can be used as a plug-in with PICProp, without affecting the propagation of aleatoric uncertainty.\n\nWe would like to clarify that we do not neglect the significance of epistemic uncertainty. Our primary goal has been to introduce a completely new problem and perspective in UQ that is accompanied by theoretical guarantees. \nThis goal has required formulating a new problem, proving theoretical guarantees, comparing with alternative approaches at a theoretical level (see Appendix C), and developing efficient methods for\nsolving the resulting bi-level optimization problem. \nWe believe that transforming our idea into a general and versatile UQ framework is a process that may require several steps/papers in the near future.\n\n**To Weakness 2:** We thank the reviewer for pointing out this crucial issue. We conducted the following experiment to address the reviewer's concern about the validity of CIs.\n\nFor the pedagogical example, we independently sampled 500 boundary conditions and propagate the corresponding Chi-squared CIs using EffiPICProp to construct CIs for $x \\in [-1, 1]$. \nAccording to the Definition of joint CIs (see Definition 3.1), an instance of CI is valid if it contains the solution $u(x)$ for any $x \\in [-1, 1]$.\nIn this regard, we uniformly sampled 100 data points from $[-1, 1]$ and regard an instance of CI as valid if the solution was enclosed within the interval for these 100 data points.\nFurther, we run vanilla PINNs with clean boundary data and use the max absolute error among these 100 data points as $\\eta$.\n\nAs a result, 486 out of these 500 (0.972 $>$ 0.95) instances of CI are valid, which indicates that our approach to constructing CI is empirically valid and conservative.\n\nDue to the limited time, we will add more results in a later revision.\n\n**To Weakness 3:** We appreciate the valuable suggestions from the reviewers regarding the formatting of the article. Due to space limitations, we have placed Algorithms 1 and 2 in the appendix. However, when reformatting, we will consider moving Algorithms 1 and 2 to the main body of the article.\n\n**To Q1:** The reviewer is right. We will correct this typo in the revised manuscript.\n\n**To Q2:** We assume the reviewer meant \"Four lines above eqn. (5)\". The reviewer is right. We will correct this typo in the revised manuscript.\n\n**To Q3:** Thanks for the reviewer's suggestion. We use the term \"SimPICProp\" for the naive regression approach to distinguish it from the meta-learning approach, EffiPICProp ($\\lambda > 0$). We will consider the reviewer's suggestion for revision.\n\n**To Q4:** By CI in Def. 3.1, we meant **joint** confidence interval, which is exactly confidence sequence, to the best of the authors' knowledge. \nWe use the term (joint) confidence interval as it is commonly used in the literature ([27][29][19]).\nFor revision, we will cite the work of Ramdas and point out the equivalence between confidence sequence and joint confidence interval. \nWe will also modify Def. 3.1 to avoid ambiguity."}]

Author rebuttal:

Meta-review: This paper addresses the important problem of uncertainty estimation and propagation in physically informed ML systems. The authors proposed a confidence interval propagation method for PDE systems. Specifically,  this paper proposes a method called Physically Informed Confidence Interval Propagation (PICProp), which propagates a confidence interval from any boundary point in a domain to any other point in the domain, and conducts a theoretical analysis. They further show that PICProp is computationally expensive and may not be applicable outside of small regions, and to address this, the authors develop an efficient version of PICProp known as EffiPICProp that makes confidence interval propagation more computationally feasible. The authors validate their method through computational experiments.
Decision: Accept (poster)