remarkRemark hypothesisHypothesis claimClaim NeuralUQZ. Zou, X. Meng, A. F.
Psaros, G. E. Karniadakis

# NeuralUQ: A comprehensive library for uncertainty quantification in neural
differential equations and operators

Zongren Zou Division of Applied Mathematics, Brown University, Providence, RI
02912, USA. Xuhui Meng Institute of Interdisciplinary Research for Mathematics
and Applied Science, School of Mathematics and Statistics, Huazhong University
of Science and Technology, Wuhan 430074, China. (Corresponding author: ).
<EMAIL_ADDRESS>Apostolos F Psaros11footnotemark: 1 George Em
Karniadakis∗ Pacific Northwest National Laboratory, Richland, WA 99354, USA.

###### Abstract

Uncertainty quantification (UQ) in machine learning is currently drawing
increasing research interest, driven by the rapid deployment of deep neural
networks across different fields, such as computer vision, natural language
processing, and the need for reliable tools in risk-sensitive applications.
Recently, various machine learning models have also been developed to tackle
problems in the field of scientific computing with applications to
computational science and engineering (CSE). Physics-informed neural networks
and deep operator networks are two such models for solving partial
differential equations and learning operator mappings, respectively. In this
regard, a comprehensive study of UQ methods tailored specifically for
scientific machine learning (SciML) models has been provided in [45].
Nevertheless, and despite their theoretical merit, implementations of these
methods are not straightforward, especially in large-scale CSE applications,
hindering their broad adoption in both research and industry settings. In this
paper, we present an open-source Python library (github.com/Crunch-UQ4MI),
termed _NeuralUQ_ and accompanied by an educational tutorial, for employing UQ
methods for SciML in a convenient and structured manner. The library, designed
for both educational and research purposes, supports multiple modern UQ
methods and SciML models. It is based on a succinct workflow and facilitates
flexible employment and easy extensions by the users. We first present a
tutorial of NeuralUQ and subsequently demonstrate its applicability and
efficiency in four diverse examples, involving dynamical systems and high-
dimensional parametric and time-dependent PDEs.

###### keywords:

education software, uncertainty quantification, scientific machine learning,
PINNs, DeepONet, deep learning

65-01, 65-04, 65L99, 65M99, 65N99

##  1 Introduction

Physics-informed and more generally scientific machine learning (SciML)
provide potent tools for accurately modeling and predicting the dynamic
response of physical systems; see [26, 13, 66] for a comprehensive review as
well as [7, 45]. Neural networks (NNs) for solving ordinary/partial
differential equations (ODEs, PDEs), such as the physics-informed NN (PINN),
and for learning operator mappings, such as the deep operator network
(DeepONet), seamlessly combine scattered observational data with available
physical models (e.g. [38, 27, 61, 52, 57]). Specifically, they can address
mixed problems when only partial observational data and information regarding
the underlying physical system are available (e.g., [51, 8]); it is
straightforward to incorporate noisy and multi-fidelity data (e.g., [70, 41,
40, 39]); they can provide physically consistent predictions even for
extrapolatory tasks (e.g., [70]); they do not require computationally
expensive mesh generation (e.g., [50, 37]); and open-source software is
available and continuously extended (e.g., [37, 23, 29, 10, 47, 48, 21, 69]).

However, uncertainty, associated with noisy and limited data as well as NN
overparameterization, for example, can degrade the accuracy of these models
significantly [2, 45, 60, 19]. As a result, uncertainty quantification (UQ) is
paramount for deep learning to be reliably used in critical applications
involving physical and biological systems [67, 46, 59, 3, 77, 62]. In Fig. 1,
we provide an overview of SciML problems of interest as well as different
sources of uncertainty. In addition to reliability concerns, UQ is also
necessary for efficient and economical design; for capital budgeting and risk-
adaptive decision making in related projects; and even for assessing
investment opportunities involving such systems in the form of real options
[53, 15]. Clearly, this broad view of applications gives to UQ an additional
dimension pertaining to the potential of harnessing uncertainty, and renders
it an indispensable research area of SciML applied to computational science
and engineering problems and beyond.

Figure 1:  Schematic overview of a SciML problem scenario. As shown in (a),
SciML problems of interest include, but are not limited to, forward, inverse
and mixed ODEs/PDEs, and operator learning problems. Next, in the pie chart of
(b) a qualitative breakdown of total uncertainty is shown, describing the
contributions from data (noisy, gappy); physical models (misspecification,
stochasticity); neural networks (architecture, hyperparameters,
overparametrization); and posterior inference. Note that we mainly address two
types of uncertainty in the current paper, i.e., aleatoric (data noise) and
epistemic (NN overparametrization and limited data) uncertainties. As shown in
(c), active learning for improving the predictions can subsequently be
performed, i.e., to acquire new inputs according to an acquisition function
that is based on the output uncertainty from part (b).

Further, UQ in SciML is considerably more challenging than in traditional
numerical solvers because of optimization-related issues of over-parameterized
NNs, as well as the various uncertainties shown in Fig. 1. In addition, SciML
is applied to increasingly more diverse problems, from steady to time-
dependent and high-dimensional ones, and to problems involving systems given
only in the form of black-box models (see, e.g., [70, 41, 45, 22, 14]).
Moreover, evaluating the performance of the available UQ methods is not
straightforward, especially in situations involving limited data. In this
regard, Psaros et al. [45] provided an extensive review and proposed novel
methods of UQ for SciML. Nevertheless, employing current methods in large-
scale applications and developing new ones, where trustworthy implementations
and fast experimentation are needed, requires an up-to-date, well-maintained,
and versatile software package. To this end, we present in this paper an open-
source Python library termed NeuralUQ, which can be accessed at
github.com/Crunch-UQ4MI. We expect that NeuralUQ will be a useful toolbox for
both practitioners and researchers, who either choose to apply off-the-shelf
methods for solving their problems, or to build their own methods based on the
ones the library currently supports. In passing, we note that NeuralUQ is
constantly updated and its most recent capabilities can be viewed in the
GitHub repository.

We organize the paper as follows. In Section 2, we briefly present SciML,
specifically the PINN and DeepONet methods, as well as recent advances in UQ
for SciML, which are the building blocks for the library; however, other
neural models for PDEs and operators can be employed as well. In Section 3, we
describe the design of the NeuralUQ library, which includes model
construction, posterior inference, and performance evaluation as well as post-
training calibration. In Section 4, we provide four computational examples for
demonstrating the applicability and capabilities of our library. Lastly, in
Section 5, we conclude with the findings of this work and discuss future
research.

##  2 Uncertainty quantification for scientific machine learning

In this section, we first briefly describe various SciML problem scenarios and
the recently developed deep learning methods for solving them. Subsequently,
we review the related UQ methods as well as the approaches for model
performance evaluation and for post-training calibration. We also provide a
glossary table of frequently used terms in this paper, and in ML/SciML more
generally, in Appendix 6.1. The interested reader is directed to [45] for more
detailed information regarding UQ methods for SciML.

###  2.1 Scientific machine learning

Consider a nonlinear ODE/PDE describing the dynamics of a physical system as
follows:

(1a) $\displaystyle\pazocal{F}_{\lambda,\xi}[u(x;\xi)]$
$\displaystyle=f(x;\xi)\text{, }x\in\Omega\text{, }\xi\in\Xi,$ (1b)
$\displaystyle\pazocal{B}_{\lambda,\xi}[u(x;\xi)]$
$\displaystyle=b(x;\xi)\text{, }x\in\Gamma,$

where $x$ is the $D_{x}$-dimensional space-time coordinate, $\Omega$ is a
bounded domain with boundary $\Gamma$, $\xi$ is a random event in a
probability space $\Xi$, while $f(x;\xi)$ and $u(x;\xi)$ represent the
$D_{u}$-dimensional source term and sought solution evaluated at $(x,\xi)$,
respectively. Further, $\pazocal{F}_{\lambda,\xi}$ is a general differential
operator; $\pazocal{B}_{\lambda,\xi}$ is a boundary/initial condition operator
acting on the domain boundary $\Gamma$; and $b$ is a boundary/initial
condition term. Furthermore, $\lambda(x;\xi)$, with
$x\in\Omega_{\lambda}\subset\Omega$, denotes the problem parameters. We refer
to problems where the operators $\pazocal{F}_{\lambda,\xi}$ and
$\pazocal{B}_{\lambda,\xi}$ are known, i.e., the differential operators that
define them are specified, as neural ODEs/PDEs, because we use neural networks
(NNs) throughout the present work to solve them. We refer to problems
involving unknown operators $\pazocal{F}_{\lambda,\xi}$ and
$\pazocal{B}_{\lambda,\xi}$ as neural operators, because we use NNs to learn
them.

For neural ODEs/PDEs, $\pazocal{F}_{\lambda,\xi}$, $\pazocal{B}_{\lambda,\xi}$
are known and deterministic, i.e. $\xi$ is fixed, while $\lambda(x)$, and
$u(x)$ are partially unknown. Given noisy data of $u$ and $\lambda$ sampled at
random finite locations in $\Omega$ and $\Omega_{\lambda}$, respectively, as
well as noisy data of $f$ and $b$, the objective of this problem is to obtain
the solution $u(x)$ at every $x\in\Omega$ and the problem parameters
$\lambda(x)$ at every $x\in\Omega_{\lambda}$. The dataset in this case is
expressed as
$\pazocal{D}=\\{\pazocal{D}_{f},\pazocal{D}_{b},\pazocal{D}_{u},\pazocal{D}_{\lambda}\\}$,
where $\pazocal{D}_{u}=\\{x_{i},u_{i}\\}_{i=1}^{N_{u}}$,
$\pazocal{D}_{\lambda}=\\{x_{i},\lambda_{i}\\}_{i=1}^{N_{\lambda}}$,
$u_{i}=u(x_{i})$, $\lambda_{i}=\lambda(x_{i})$, and
$\pazocal{D}_{f},\pazocal{D}_{b}$ are given as above. The PINN method,
developed in [50], addresses the PDE problems deterministically, by
constructing NN approximators, parameterized by $\theta$, for modeling the
solution $u(x)$. The approximator $u_{\theta}(x)$ is substituted into Eq. (1)
via automatic differentiation [1] for producing
$f_{\theta}=\pazocal{F}_{\lambda,\xi}[u_{\theta}]$ and
$b_{\theta}=\pazocal{B}_{\lambda,\xi}[u_{\theta}]$. For solving the mixed
deterministic PDE problem, an additional NN can be constructed for modeling
$\lambda(x)$. Finally, the optimization problem for obtaining $\theta$ in the
mixed problem scenario with the dataset
$\pazocal{D}=\\{\pazocal{D}_{f},\pazocal{D}_{b},\pazocal{D}_{u},\pazocal{D}_{\lambda}\\}$
is cast as

(2)
$\hat{\theta}=\underset{\theta}{\mathrm{argmin}}~{}\pazocal{L}(\theta)\text{,
where
}\pazocal{L}(\theta)=\frac{w_{f}}{N_{f}}\sum_{i=1}^{N_{f}}||f_{\theta}(x_{i})-f_{i}||_{2}^{2}\\\
+\frac{w_{b}}{N_{b}}\sum_{i=1}^{N_{b}}||b_{\theta}(x_{i})-b_{i}||_{2}^{2}+\frac{w_{u}}{N_{u}}\sum_{i=1}^{N_{u}}||u_{\theta}(x_{i})-u_{i}||_{2}^{2}+\frac{w_{\lambda}}{N_{\lambda}}\sum_{i=1}^{N_{\lambda}}||\lambda_{\theta}(x_{i})-\lambda_{i}||_{2}^{2}.$

In this equation, $\\{w_{f},w_{b},w_{u},w_{\lambda}\\}$ are objective function
weights for balancing the various terms in Eq. (2). Following determination of
$\hat{\theta}$, $u_{\hat{\theta}}$ and $\lambda_{\hat{\theta}}$ can be
evaluated at any $x$ in $\Omega$, $\Omega_{\lambda}$, respectively.

For neural operators, $\pazocal{F}_{\lambda,\xi}$, $\pazocal{B}_{\lambda,\xi}$
are unknown and stochastic, while $\lambda(x;\xi),u(x;\xi)$ are partially
unknown. Note that the stochasticity of $\pazocal{F}_{\lambda,\xi}$ and
$\pazocal{B}_{\lambda,\xi}$ arises from the stochasticity in $\lambda(x;\xi)$
with $\xi\in\Xi$. The available data consists of paired realizations of
$f,b,u$, and $\lambda$, each pair associated with a different $\xi\in\Xi$. The
realizations are considered noisy in general for $u$, and clean for $f,b$, and
$\lambda$. Further, the objective of this problem is to learn the operator
mapping defined by Eq. (1). The DeepONet method, developed in [35] and
designed for general nonlinear operator learning, addresses this problem by
constructing an approximator using NNs, denoted by $G_{\theta}$ and
parameterized by $\theta$, to model the mapping from the input function domain
to the output function domain. In the following and without loss of
generality, we consider both $f$ and $b$ in Eq. (1) as known and the problem
becomes learning the mapping from the parameter $\lambda$ to the solution $u$.
Further, each value of $\xi\in\Xi$ in Eq. (1) corresponds to a different
$\lambda(x;\xi)$ function, and thus to a different solution $u(x;\xi)$. For
any $\xi\in\Xi$, the approximator $G_{\theta}$ takes as input a representation
of $\lambda(x;\xi),x\in\Omega_{\lambda}$, and outputs a function,
$G_{\theta}(\lambda)(x;\xi),x\in\Omega$. The DeepONet consists of two sub-
networks, one that takes the coordinates $x\in\Omega$ as input, called trunk
net, and one that takes the discretized function $\lambda$ as input, called
branch net. Considering a one-dimensional $u$, i.e., $D_{u}=1$, each of the
two networks produces a $w$-dimensional output. The two outputs are merged by
computing their inner product for producing the final prediction
$G_{\theta}(\lambda)(x),x\in\Omega$. For learning the mapping from $\lambda$
to $u$ using the vanilla DeepONet [36], the minimization problem with the mean
squared error (MSE) loss

(3)
$\hat{\theta}=\underset{\theta}{\mathrm{argmin}}~{}\pazocal{L}(\theta)\text{,
where
}\pazocal{L}(\theta)=\frac{1}{NN_{u}}\sum_{i=1}^{N}\sum_{j=1}^{N_{u}}||G_{\theta}(\lambda(.;\xi_{i}))(x_{j})-u_{j}^{(i)}||_{2}^{2}$

is solved, using the training dataset
$\pazocal{D}^{\prime}=\\{x_{j},\xi_{i},\lambda^{(i)},u_{j}^{(i)}\\}_{i,j=1,1}^{N,N_{u}}$,
where $\\{x_{j}\\}_{j=1}^{N_{u}}$ are the locations, $\\{\xi_{i}\\}_{i=1}^{N}$
are the sampled events and in this case represent the unique pairing between
$u$ and $\lambda$, $u_{j}^{(i)}:=u(x_{j},\xi_{i})$ is the measurement of $u$
at the location $x_{j}$ for event $\xi_{i}$, and $\lambda^{(i)}$ is the
representation of $\lambda(x;\xi_{i})$. Note that, to be consistent with
neural ODEs/PDEs, we use $u_{\theta}$ to denote the approximator for operators
in the rest of this paper. In this case,
$u_{\theta}=u_{\theta}(x;\lambda(\cdot;\xi))$. DeepONets, once trained, could
act as PDE-agnostic surrogates encoding the hidden physics for further
learning and inferences [41]. Finally, note that in addition to supervised
learning using data as described above, DeepONet can alternatively be trained
by employing only the PINN loss during training. As shown in [20], the most
effective training is the hybrid one, where both physics and data are used in
the training loss; see also [64].

###  2.2 Uncertainty quantification

In contrast to the point estimates obtained via the optimization problems of
Eqs. (2) and (3), the goal of UQ is to estimate the predictive distribution

(4) $\displaystyle p(u|x,\pazocal{D})=\int
p(u|x,\theta)p(\theta|\pazocal{D})d\theta,$

which is often referred to as Bayesian model average (BMA). In Eq. (4),
$p(\theta|\pazocal{D})$ denotes the posterior density, which can be obtained
using Bayes’ rule via

(5)
$p(\theta|\pazocal{D})=\frac{p(\pazocal{D}|\theta)p(\theta)}{p(\pazocal{D})}.$

In Eq. (5), $p(\pazocal{D}|\theta)$ is the likelihood of the data, i.e.,
$p(\pazocal{D}|\theta)=\prod_{i=1}^{N}p(u_{i}|x_{i},\theta)$ for independent
and identically distributed (i.i.d.) data; $p(\theta)$ is the probability
density function for the prior distribution of the parameters $\theta$; and
$p(\pazocal{D})$ is called marginal likelihood or evidence, because it
represents the probability that we observe $\pazocal{D}$. The evidence is
given as

(6) $p(\pazocal{D})=\int p(\pazocal{D}|\theta)p(\theta)d\theta.$

The uncertainty of parameters $\theta$ is typically referred to as epistemic
uncertainty.

Figure 2: A unified view of the UQ models studied in Section 2.2; adopted from
[45]. The model $\pazocal{M}_{1}$ can be a NN or a generator network, or
nothing. The model $\pazocal{M}_{2}$ can be a PDE or an operator, such as a
pre-trained DeepONet. The quantities $\phi$, $\psi$ can be the solution $u$
and the source $f$ in Eq. (1), or vice-versa; $x$ represents the space-time
input domain; and $\theta$ represents the parameters of the NN or the input to
the generator. Further, “$\chi$” is used for modeling additional fields and
can represent, for example, the model parameters $\lambda$ in Eq. (1). The
dashed arrow indicates that the connection may not exist for some problem
scenarios.

Generally, inferring $\theta$ given $\pazocal{D}$, i.e., computing the
posterior exactly via Eq. (5) is computationally and analytically intractable.
This is addressed by approximate posterior inference, which approximates the
posterior by another distribution and/or obtains representative samples from
the posterior. Specifically, the Bayesian and the deterministic frameworks are
the two main families of UQ methods for SciML. The former employs either a
sampling method or a variational inference (VI) method to obtain samples from
the posterior distribution, whereas the latter can be interpreted as an
extension of standard training of deterministic NNs. Various UQ methods have
been developed for SciML within these two frameworks. We adopt the unified
framework proposed in [45] to include all the UQ models employed in this study
in Fig. 2, and present a short overview for them in the rest of this section.

For neural ODEs/PDEs, the solution $u$ and/or the problem parameters $\lambda$
are modeled with surrogates parameterized by $\theta$ and denoted by
$\pazocal{M}_{1}$ in Fig. 2. The Bayesian PINNs (B-PINNs) [72] and the
functional prior (FP) based on the physics-informed generative adversarial
networks (PI-GANs) [41, 71] are two recently proposed physics-informed
learning approaches with UQ within the Bayesian framework. The Bayesian neural
networks (BNNs) and the generators from trained PI-GANs are employed as the
surrogate for $u/\lambda$ in these two methods, respectively. Further, either
sampling or variational inference (VI) approaches can be utilized for
obtaining the posterior samples for $\theta$. In addition to the
aforementioned methods, the Monte Carlo dropout [18] (MCD) as well as the deep
ensemble [31] (DEns) approaches, which employ standard, i.e., deterministic,
NNs as surrogates for $u$ and/or $\lambda$, can be combined with PINNs to
quantify uncertainties in problems with known physics [76, 45].

For neural operators, a DeepONet parameterized by $\theta$ is employed as the
surrogate for the solution $u$ and denoted by $\pazocal{M}_{2}$ in Fig. 2. In
the Bayesian framework, [34] employed BNNs in DeepONet and a sampling method
is then utilized to obtain the posterior samples for the parameters in BNNs.
As for deterministic methods, extensions of standard training of DeepONets,
such as DEns, are used for quantifying uncertainties in predictions [45, 73,
43]. These models correspond to Fig. 2 with $\pazocal{M}_{2}$ being a DeepONet
and the dashed arrows are not used. Further, the inference input in these
models, i.e. the input to the branch net, should be represented in the same
way as used in the training, i.e., the number of the discrete values as well
as their locations should be the same. For a more general case in which the
inference input and/or output data are noisy and measured at a certain number
of random locations, [41] proposed to combine a pre-trained DeepONet with a
GAN-based FP to quantify the uncertainties in predictions. In [45], a BNN is
utilized as an alternative to the generator of a GAN in [41]. These two models
correspond to Fig. 2 with $\pazocal{M}_{1}$ being either a BNN or the
generator from a pre-trained GAN and $\pazocal{M}_{2}$ a pre-trained DeepONet.

In all cases, different $\theta$ values correspond to different functions
$u_{\theta}$. We refer to the surrogate combined with a distribution
$p(\theta)$ for $\theta$ as a process. Following approximate posterior
inference, the posterior samples for $\theta$ are used for obtaining the
statistics of the function samples. In this regard, the value of $u$ at a
location $x$ given the data $\pazocal{D}$ is a random variable denoted by
$(u|x,\pazocal{D})$. The mean of $(u|x,\pazocal{D})$ is denoted by
$\hat{u}(x)=\mathbb{E}[u|x,\pazocal{D}]$ and approximated by $\bar{\mu}(x)$ as

(7)
$\hat{u}(x)\approx\bar{\mu}(x)=\frac{1}{M}\sum_{j=1}^{M}u_{\hat{\theta}_{j}}(x),$

where $\\{u_{\hat{\theta}_{j}}(x)\\}_{j=1}^{M}$ is the set of NN predictions
corresponding to the samples $\\{\hat{\theta}_{j}\\}_{j=1}^{M}$. Further, the
diagonal part of the covariance matrix of $(u|x,\pazocal{D})$ is given by

(8)
$Var(u|x,\pazocal{D})\approx\bar{\sigma}^{2}(x)=\underbrace{\Sigma_{u}^{2}}_{\bar{\sigma}_{a}^{2}(x)}+\underbrace{\frac{1}{M}\sum_{j=1}^{M}(u_{\hat{\theta}_{j}}(x)-\bar{\mu}(x))^{2}}_{\bar{\sigma}_{e}^{2}(x)}.$

In this equation, $\bar{\sigma}^{2}(x)$ represents the approximate total
uncertainty of $(u|x,\pazocal{D})$, while $\bar{\sigma}_{a}^{2}(x)$ and
$\bar{\sigma}_{e}^{2}(x)$ denote the aleatoric and the epistemic parts of the
total uncertainty, respectively. In practice, a Gaussian (or, e.g., a
Student-t) distribution is commonly fitted to the posterior samples. For
example, we denote as $\pazocal{N}(\bar{\mu}(x),\bar{\sigma}^{2}(x))$ the
Gaussian approximation to the samples
$\\{u_{\hat{\theta}_{j}}(x)\\}_{j=1}^{M}$ at an arbitrary location $x$, where
$\bar{\mu}(x)$ and $\bar{\sigma}^{2}(x)$ are obtained via Eqs. (7) and (8),
respectively.

###  2.3 Performance evaluation and post-training calibration

Performance evaluation is necessary for model selection, e.g.,
prior/architecture optimization, comparison between UQ methods, and overall
quality evaluation of the UQ design procedure [44]. As reviewed in [45],
performance evaluation metrics include, but are not limited to, the relative
$\ell_{2}$ error for evaluating accuracy, the mean predictive likelihood for
evaluating predictive capacity, and the root mean squared calibration error
for evaluating statistical consistency; i.e., how closely the predicted model
matches the data-generating distribution.

Further, calibration approaches can be adopted to improve statistical
consistency in UQ for SciML [45]. Various approaches have been developed for
calibration, e.g., [33] re-weights the predicted variance optimally to
calibrate the output uncertainty; [30] approximates the mis-calibration
function and applies it to all predictive distributions, and thus modifies
both the output mean and variance; and [75] fits an empirical distribution to
the scaled residuals, and also calibrates both the output mean and variance.
To use these approaches, an additional dataset is required after training.
This can be a left-out calibration dataset, similar to validation and test
sets.

##  3 Overview and details of NeuralUQ

In this section, we first present an overview of the NeuralUQ library, and
subsequently provide details regarding its design, use, and also customization
for extended capabilities.

###  3.1 Library overview

Figure 3: Library design overview. In the design phase, a model is constructed
that includes processes representing a proposed solution family and
likelihoods or loss functions representing the data distribution. Processes
include surrogates, which can be construed as assumed functional forms, and
priors distributions or variables for the parameters of the surrogates,
depending on the type of method used. Likelihoods or loss functions are
constructed using the data and the underlying neural ODE/PDE. To select the
solution function samples that best explain the data, a posterior inference
method is combined with the model. Finally, the obtained function samples can
be calibrated using additional data. In the revision phase, metrics are used
for evaluating the performance and re-designing the solution procedure.
Procedure 1 NeuralUQ for solving SciML problems.

  1. Step 1

Specify a process, which is composed of a surrogate as well as its parameters
$\theta$ to approximate the solution of an ODE/PDE, using the surrogates,
variables, and process modules.

  2. Step 2

Specify the differential equations using automatic differentiation.

  3. Step 3

Construct the likelihoods or loss functions based on the equations,
boundary/initial conditions, and/or given data, using the likelihoods module.

  4. Step 4

Combine the process and likelihoods/loss functions to construct a model.

  5. Step 5

Specify an inference method using the inference module, and compile the model
with it to obtain a sampler for posterior distributions by calling
model.compile.

  6. Step 7

Call model.run to run the compiled sampler to obtain posterior samples.

  7. Step 8

Compute predictions with uncertainties based on Eqs. (7) and (8) using the
posterior samples of $\theta$ from Step 7.

  8. Step 9

Use the metrics module for performance evaluations using the predictions from
Step 8.

  9. Step 10

Use the calibrations module for post-training calibration (if applicable).

As mentioned in the previous section, the value of $u$ at a location $x$ given
the data $\pazocal{D}$ is a random variable denoted by $(u|x,\pazocal{D})$.
The aim of UQ in the library is to obtain, calibrate, compute the statistics,
and evaluate the quality of samples drawn from the posterior distribution
$p(u|\pazocal{D})$. To that end, in the design phase a model is constructed
that includes surrogates representing the solution $u$ and potentially $f$
and/or $\lambda$, as well as likelihood functions representing the data
distribution. Surrogates are parameterized by $\theta$ and when combined with
probability distributions of $\theta$ are represented as stochastic processes;
i.e., a proposed family of solutions. To obtain the function samples for the
solutions, a posterior inference method is combined with the model. Finally,
the obtained function samples can be calibrated using additional data. In the
revision phase, metrics are used for evaluating the performance and re-
designing the solution procedure. A schematic workflow of the library is
provided in Fig. 3 and Procedure 1.

###  3.2 Library design

UQ methods in NeuralUQ
---
UQ framework | UQ family | Process | Inference methods
Bayesian | Samplable | BNNs, Gens | HMC, LD, MALA, NoUTurn
Variational | BNNs, Gens | MFVI, MCD
Deterministic | Trainable | FNNs, Gens, DeepONets | DEns, SEns, LA
Table 1:  Summary of the UQ methods implemented in NeuralUQ. NeuralUQ includes
three different UQ method families, i.e., Samplable, Variational, and
Trainable. Each UQ method family has various processes and inference methods,
e.g., BNNs+HMC in Samplable; BNNs+MFVI Variational; and FNN+DEns in Trainable.
Gens: Generators in GANs. Short and detailed descriptions for these inference
methods are provided in Table 2 and [45, 55, 68, 24, 5], respectively. Note
that BNNs/Generators/FNNs combined with known ODEs/PDEs are used for UQ for
neural ODE/PDE problems, while BNNs/Generators combined with pre-trained
DeepONets and DeepONets using DEns are used for UQ in neural operator
problems.

As discussed in Section 2.2, there are two general UQ frameworks for SciML,
namely the Bayesian and the deterministic frameworks. The former relies on
Bayes’ theorem, and includes sampling methods, such as Monte Carlo Markov
Chain (MCMC), and VI methods, such as mean-field variational inference (MFVI).
The latter relies on standard NN training, i.e., minimizing/maximizing an
objective function, typically using the MSE or the maximum a posteriori (MAP)
metric, between measurements and NN predictions.

The procedural differences of the aforementioned methods necessitate the
development of three different UQ method families in our library, i.e., the
_Samplable_ , _Variational_ , and _Trainable_ to include all the sampling-,
VI- and standard-NN-based methods, respectively. Specifically, the sampling
methods target directly the unnormalized posterior density function, without
employing and optimizing any trainable variables. The VI-based methods employ
pre-defined distributions parameterized by trainable variables, which are
optimized to approximate the true posterior distributions. Finally, the
standard-NN-based UQ methods obtain point estimates for the parameters in the
surrogates using an optimizer, rather than directly sampling from or
approximating the posterior distributions. A summary of the currently
supported UQ methods in NeuralUQ is provided in Table 1.

###  3.3 Details of NeuralUQ components

In this section, we describe the main components, built-in functionalities and
how to use NeuralUQ.

####  3.3.1 Processes

Figure 4:  A breakdown of the Process module. A Process instance contains a
surrogate as well as the attributes for its parameters. Processes, used in
different UQ method families, have different surrogates and attributes for
their parameters, which are denoted with different line colors in this figure.
Different combinations of the components yield different processes, e.g., an
fnn surrogate combined with a prescribed prior distribution for its parameters
results in a BNN.

As described in Procedure 1, the first step for using NeuralUQ is to specify a
process for the sought solution, e.g., BNNs or pre-trained generators,
according to the task at hand. In this regard, the surrogates, variables, and
Process modules have been developed in NeuralUQ. Specifically, the surrogates
module expresses how $x$ and $\theta$ are combined to yield the approximators
$u_{\theta}(x)$ or $\lambda_{\theta}(x)$ ($\pazocal{M}_{1}$ or
$\pazocal{M}_{2}$ in Fig. 2); the variables module is used to specify the
attributes of the parameters in surrogates, e.g., prior/posterior
distributions (Fig. 4); and the Process module is an interface for the
implementation of the unified models in Fig. 2, which connects the surrogate
with its associated parameters.

Currently, NeuralUQ provides four built-in surrogates, namely,
surrogates.identity, used as a surrogate for an unknown constant; fnn, for a
fully-connected neural network (FNN); generator, for the generator of a PI-GAN
pre-trained with historical data; and deeponet, for a DeepONet. The variables
module, on the other hand, is used to specify the prior and/or posterior
distributions of surrogate parameters. For instance, the prior distribution is
required for surrogate parameters in the Samplable family, while both the
prior and the posterior distributions are required for parameters in the
Variational family. The built-in variables of NeuralUQ are variables.const,
used for parameters of unknown constants as well as the generators of pre-
trained GANs/PI-GANs; and fnn, used for FNNs and also DeepONets. Further, as
presented in Fig. 4 (also in Table 1 and Section 3.2), different UQ method
families support different surrogates and are compatible with different
variables. In NeuralUQ, variables.fnn and variables.const are decomposed into
three different sub-modules, each one pertaining to a UQ method family of
Table 1, i.e., variables.fnn/const.Samplable, Variational, and Trainable.

With the specified surrogate and parameters, we can then instantiate a Process
object to create a corresponding process. The surrogate and the attributes for
its parameters can then be accessed via Process.surrogate and
Process.prior/posterior, respectively. Also, for problems with unknown
parameters or multi-physics, more than one process may be employed. The module
Process.key can be used to trace each process in NeuralUQ.

####  3.3.2 Likelihoods

In the Samplable and Variational families, likelihood distributions for the
given data are required to construct the unnormalized posterior distributions
of the surrogate parameters. Similarly, a metric, namely, MSE or MAP, is
required for methods in the Trainable family. In NeuralUQ, the likelihoods
module can be used to define the likelihood distributions in Samplable and
Variational, as well as the MSE/MAP metrics in Trainable (i.e., Step 3 in
Procedure 1). In addition, the currently supported likelihood for Samplable
and Variational is the Gaussian distribution, accessed via likelihoods.Normal,
and the supported metric for Trainable is the MSE, accessed via
likelihoods.MSE. Each likelihood receives inputs, targets and processes as
input arguments, for specifying the locations of the measurements, the
measurements, and the corresponding processes used to model them. For tasks
involving PDEs, the argument pde to the likelihood, which is a Python callable
defined separately (i.e., Step 2 in Procedure 1), is also required.

####  3.3.3 Inferences

Uncertainty in the sought solution and/or unknown problem parameters is
associated with the uncertainty of the surrogate parameters. The inferences
module has been developed to obtain the posterior samples for the surrogate
parameters in NeuralUQ (i.e., Step 5 in Procedure 1). As shown in Table 1,
NeuralUQ provides various built-in inference methods in each of the UQ method
families. Specifically, NeuralUQ currently supports inferences.HMC, LD, MALA,
and NoUTurn in Samplable; inferences.VI, MCD in Variational; and
inferences.DEns and SEns in Trainable. Short descriptions of these methods can
be found in Table 2.

Inference methods for posterior estimation
---
Full name | Abbreviation | Short description
Hamiltonian | HMC | Simulates Hamiltonian dynamics
Monte Carlo [42] | followed by the Metropolis-Hasting algorithm
Metropolis-adjusted | MALA | Simulates Langevin dynamics
Langevin algorithm [55, 68] | followed by the Metropolis-Hasting algorithm
Langevin dynamics [42, 65] | LD | One leap-frog step of HMC
NoUTurn [24, 5] | NUTS | A variant of HMC that adaptively chooses
the number of steps for the leap-frog integrator
Mean-field | MFVI |  Approximates $p(\theta|\pazocal{D})$ by a Gaussian variational
variational inference [6] | distribution with diagonal covariance matrix
Monte Carlo | MCD |  Approximates $p(\theta|\pazocal{D})$ by a Bernoulli variational
dropout [18] | distribution with fixed dropout rate
Deep ensemble [31] | DEns |  Standard training $M$ times independently
with random initializations
Snapshot ensemble [25] | SEns | Standard training with learning rate
schedule and snapshots
Laplace | LA | Standard training followed by fitting a
approximation [58, 54] |  Gaussian to approximate $p(\theta|\pazocal{D})$
Table 2:  Overview of the inference methods in NeuralUQ; adopted from [45].
HMC, MALA, LD, NUTS, and MFVI are inference methods corresponding to the
Bayesian framework in Table 1; and DEns to the deterministic one. Note that
MCD can be interpreted and implemented in both frameworks [18]. In the present
study, we treat it as a Bayesian method.

####  3.3.4 Performance evaluation and post-training calibration

Performance evaluation, model revision and post-training prediction
improvement are crucial for reliably deploying UQ methods in risk-sensitive
SciML settings. In NeuralUQ, we have developed the metrics module, which
currently includes two different metrics, i.e., metrics.RL2E/MSE, and NLL
(i.e., Step 9 in Procedure 1). The former serves as an evaluation metric of
the predicted mean and computes the relative $L_{2}$ error (RL2E) or MSE
between the predicted mean and the reference solution, whereas the latter
serves as a metric of the predicted uncertainty and computes the negative log
likelihood (NLL) [74, 49].

Further, we have developed the calibrations module for performing post-
training calibration (i.e., Step 10 in Procedure 1). Specifically, this module
currently includes one sub-module, i.e., calibrations.var, for the
implementation of the method proposed in [33]. We present a workflow for
calibration in Procedure 2, and ; a demonstration example involving post-
training calibration can be found in the GitHub repository of the library. We
will keep updating NeuralUQ to support more metrics and calibration methods,
e.g., [12], in the future.

Procedure 2 A workflow for calibration using NeuralUQ.

  1. Step 1

Obtain uncalibrated posterior function samples at the locations $x$ where new
measurements are available, using the posterior samples of the parameters from
Step 7 in Procedure 1.

  2. Step 2

Choose and apply a calibration method from the calibrations module, which
takes as input the function samples from Step 1 and new measurements.

  3. Step 3

Obtain calibrated function samples.

###  3.4 Customization

The components of NeuralUQ, including surrogates, variables, likelihoods, and
inferences, are highly flexible and customizable. In this section, we describe
how to configure the library for extending current capabilities and addressing
problem scenarios beyond the ones considered in this study.

####  3.4.1 Surrogates

As mentioned in Section 3.3.1, NeuralUQ currently supports four surrogates,
namely, identity, fnn, generator, and deeponet. Users may need to develop and
incorporate in the library alternative surrogates. In NeuralUQ, this can be
achieved in a straightforward manner by following Procedure 3. Demonstration
examples involving customizing new surrogates are presented in Appendix 6.2.

Procedure 3 Customization of a new surrogate, MySurrogate, which inherits from
surrogates.Surrogate.

⬇

class MySurrogate(Surrogate):

def __call__(self, inputs, var_list):

"""

Returns the outputs of the function, defined by samples stored in ‘var_list‘,
to the input ‘inputs‘, and also the reshaped and tiled input, for derivative
computation.

"""

####  3.4.2 Variables

As mentioned in Section 3.3.1, variables.fnn/const is decomposed into three
sub-modules for different UQ method families, and can be customized according
to the task at hand. In this section, customization of the variables in the
first two families is discussed, as the corresponding variables influence more
the performance as compared to the variables in the third one. The built-in
prior distribution for Samplable and Variational variables is the Gaussian
distribution. A straightforward customization in this context pertains to
modifying the mean and/or the standard deviation to obtain different priors.
Indicatively, adding a Samplable in NeuralUQ variable can be done by following
Procedure 4, which includes the setup of the initial values as well as the log
probability density function. Another option pertains to selecting alternative
parameterized distributions for performing VI. Demonstration examples
involving prior customization as well as selecting alternative distributions
for performing VI are presented in Appendix 6.3.

Procedure 4 Customization of a new Samplable variable, MySamplable, which
inherits from variables._Samplable.

⬇

class MySamplable(_Samplable):

@property

def initial_values(self):

"""Returns initial values of the variable."""

def log_prob(self, samples):

"""Returns the log probability density of ‘samples‘."""

####  3.4.3 Inference methods

As presented in Section 3.3.3 and summarized in Table 2, NeuralUQ supports
various inference methods for each of the UQ method families, i.e., Samplable,
Variational, and Trainable. However, depending on the task at hand, users may
need to develop in the library new inference methods, e.g., for reducing
computational cost. Indicatively, adding a new sampling method can be done by
following Procedure 5. In passing, note that all sampling methods included
currently in NeuralUQ are implemented using the TensorFlow Probability MCMC
python package [32, 16], which is computationally efficient and
straightforward to use. For consistency reasons, we recommend that users of
NeuralUQ employ this package for adding new sampling methods to the library.

Procedure 5 Customization of a new inference, MyInference, which inherits from
inferences.Inference.

⬇

class MyInference(Inference):

def make_sampler(self):

"""Creates and returns a sampler, which will be used in ‘sampling‘ method."""

def sampling(self):

"""Performs the inference and returns samples."""

##  4 Demonstration examples

In this section, we provide four computational examples for demonstrating the
applicability and capabilities of our library. Specifically, we quantify
uncertainty in both neural ODEs/PDEs and neural operators applied to the
following problems: (1) a Kraichnan-Orszag system, (2) a susceptible-infected-
recovered-dead (SIRD) model for COVID-19, (3) a Korteweg-de Vries problem, and
(4) a 100-dimensional Darcy problem. In the first three cases, we obtain and
present results using uncertain PINNs, i.e., B-PINNs combined with HMC and
PINNs combined with DEns, as representative examples of the Bayesian and
deterministic methods, respectively. For the 100-dimensional Darcy problem, we
present results using the DeepONet combined with two physics-agnostic models,
i.e., physics-agnostic GAN functional prior (PA-GAN-FP) and physics-agnostic
BNN functional prior (PA-BNN-FP), as well as DEns. Results corresponding to
additional UQ methods and details on the computations in each case, e.g.,
architectures of NNs and parameters used in the training, are provided in the
GitHub repository of our library.

###  4.1 Kraichnan-Orszag system

The Kraichnan-Orszag model consisting of three coupled nonlinear ODEs
describes the temporal evolution of a system with several interacting inviscid
shear waves, whose initial conditions are drawn from a Gaussian distribution
[63]. Here, we use NeuralUQ to solve the inverse Kraichnan-Orszag problem,
i.e., identify the unknowns in the governing equations given sparse and noisy
measurements on the solutions without the knowledge of the initial conditions.
The governing equations are:

(9) $\displaystyle\frac{dx_{1}}{dt}$
$\displaystyle=ax_{2}x_{3},~{}\frac{dx_{2}}{dt}=bx_{1}x_{3},~{}\frac{dx_{3}}{dt}=-(a+b)x_{1}x_{2},$
(10) $\displaystyle x_{1}(0)$ $\displaystyle=x_{2}(0)=1,~{}x_{3}(0)=0.5,$

where $a=1$ and $b=1$ are constants.

Figure 5:  Inverse Kraichnan-Orszag problem: Predictions corresponding to
$x_{1}$, $x_{2}$, and $x_{3}$ using HMC (a) and DEns (b). Epistemic
uncertainty is denoted by $2\bar{\sigma}_{e}$, whereas total uncertainty by
$2\bar{\sigma}$.

Particularly for this problem, we assume that $a$ and $b$ are unknown, whereas
11 noisy measurements of $x_{1}$ and $x_{3}$ are available, and 7 of $x_{2}$;
see also Fig. 5. The noise for all measurements is assumed to be Gaussian with
zero mean and a known standard deviation set to 0.05. The objective here is to
identify $a$ and $b$, as well as to reconstruct the complete solutions
$x_{1}$, $x_{2}$, and $x_{3}$ for $t\in[0,10]$, with uncertainties.

We present the predictions for $x_{1},x_{2},x_{3}$ and $a,b$ using different
approaches in Fig. 5 and Table 3, respectively. As shown in Fig. 5, for both
HMC and DEns, the predicted uncertainty increases at locations with no
available data. Further, the error between the predicted mean and the
reference solution for $x_{1},x_{2},x_{3}$ is bounded by the total
uncertainty. Furthermore, as shown in Table 3, the predicted values for $a$
and $b$, using both approaches, are close to the reference values. In the
table, the results using MFVI and MCD are also provided for comparisons. In
Table 4, we evaluate the predictions for $x_{1},x_{2},x_{3}$ using the RL2E
and the NLL metrics. As shown, HMC and DEns provide the most reasonable
predictions for the mean as well as uncertainties according to the computed
metrics.

| HMC | DEns | MFVI | MCD
---|---|---|---|---
$a$ (mean $\pm$ std) | $0.8380\pm 0.0813$ | $0.8720\pm 0.0812$ | $0.5706\pm 0.0205$ | $0.6141\pm 0.0202$
$b$ (mean $\pm$ std) | $1.0705\pm 0.0454$ | $1.1463\pm 0.0558$ | $0.9388\pm 0.0166$ | $0.8931\pm 0.0161$
Table 3:  Inverse Kraichnan-Orszag problem: Predictions corresponding to $a$ and $b$ using different methods. The reference values are $a=b=1$.  | HMC | DEns | MFVI | MCD
---|---|---|---|---
RL2E | 0.0579 | 0.0585 | 0.2232 | 0.1527
NLL ($\times 10^{2}$) | -3.6558 | -4.0275 | -0.1187 | -1.8930
Table 4:  Inverse Kraichnan-Orszag problem: performance evaluation of the
considered methods on test dataset. Note that NLL is computed using the total
uncertainty. Smaller RL2E and NLL mean better accuracy for the predictions of
mean and uncertainty, respectively. RL2E corresponds to relative $L_{2}$
error, whereas NLL to the negative log likelihood.

###  4.2 SIRD model for COVID-19

Scientific models are critical tools for understanding and forecasting the
spread of the COVID-19. However, constructing models that are capable of
accurately describing the COVID-19 still remains a challenge. Here, we
demonstrate the capability of NeuralUQ for learning the corresponding model
with uncertainties given measurements. In particular, we consider the SIRD
compartmental model with time-dependent parameters/coefficients/rates for the
COVID-19 [11, 9, 56], expressed as:

(11) $\displaystyle\frac{dS}{dt}$
$\displaystyle=\frac{-\beta(t)SI}{n},~{}\frac{dI}{dt}=\frac{\beta(t)SI}{n}-\gamma(t)I-\nu(t)I,$
$\displaystyle\frac{dR}{dt}$
$\displaystyle=\gamma(t)I,~{}\frac{dD}{dt}=\nu(t)I,$

where $\beta$, the transmission rate, $\gamma$, the recovery rate, and $\nu$,
the death rate, are unknown and assumed to be time-dependent functions, and
$n$ is the total population, which is assumed to be constant and equal to the
summation of $S,I,R,D$ throughout the time span. $S(t),I(t),R(t)$ and $D(t)$
are the number of individuals susceptible to the disease, the number of
individuals infected, the (cumulative) number of individuals recovered, and
the (cumulative) number of individuals deceased from the disease,
respectively, at time $t$. In this example, we use the COVID-19 dataset from
Italy [17], which is the same as in [28].

Figure 6:  Inverse SIRD problem: Predictions for $I$, $R$, $D$ and also
$\beta$, $\gamma$ and $\nu$, using HMC (a) and DEns (b). Total uncertainty is
denoted by $2\bar{\sigma}$.

In this problem, we utilize the data on $I$, $R$, and $D$ in the first 160
days for training. The objective is to predict $I$, $R$, and $D$ for
$t\in[161,200]$, and identify the unknown functions $\beta$, $\gamma$, and
$\nu$ for $t\in[0,200]$. The results from HMC and DEns are presented in Fig.
6. We observe that the error between the predicted mean and the reference
solution of $I/R/D$ from both approaches is bounded by the total
uncertainties. Further, HMC and DEns yield different predictions for $\beta$
and $\nu$ for times $t<100$, while they are quite similar for $t>100$. Also,
both methods predict a decrease of $\beta$ and $\nu$ with time. In addition,
the predictions for $\gamma$ from both methods are similar, i.e., they are
close to zero for $t\in[0,200]$. Note that reference solutions for the unknown
functions $\beta$, $\gamma$, and $\nu$, are not available and thus
corresponding comparisons are not provided here.

Figure 7:  Inverse KdV problem: Predictions for $u(x,t)$ at different times,
as obtained by using HMC and DEns. (a) Locations of training data for $u$ and
$f$; black circles are used for $u$ and magenta cross symbols for $f$. (b)
HMC; $t=-1.5$ in the left side and $t=0.4$ in the right side. (c) DEns;
$t=-1.5$ in the left side and $t=0.4$ in the right side. Epistemic uncertainty
is denoted by $2\bar{\sigma}_{e}$, whereas total uncertainty by
$2\bar{\sigma}$. Note that the scaling for $t$ is the same as in [4]. One
could also use a different scaling for $t$ to make it positive only.

###  4.3 Inverse KdV problem

The Korteweg-de Vries (KdV) equation, derived more than 150 years ago, can be
used to model surface waves in a shallow canal. It has been applied in various
disciplines for modeling surface gravity waves, internal solitons in the
ocean, magma flows, and conduit waves, for example. In this section, we solve
an inverse KdV problem using NeuralUQ. Specifically, we discover the KdV
equation given data for a system with interactions between two solitary waves.
The KdV equation is expressed in the general form:

(12) $u_{t}-\lambda_{1}uu_{x}-\lambda_{2}u_{xxx}=f,$

where $\lambda_{1}$ and $\lambda_{2}$ are constants, and $f=0$ is considered.
The exact solution can be expressed as [4]

(13) $u=2\partial^{2}_{x}R,$

where

(14)
$R=\log\left[\exp(-\eta_{1}-\eta_{2})+\exp(\eta_{1}-\eta_{2})+\exp(\eta_{2}-\eta_{1})+\frac{(a_{1}-a_{2})^{2}}{(a_{1}-a_{2})^{2}}\exp(\eta_{1}+\eta_{2})\right],$

and

(15) $\displaystyle\eta_{1}$ $\displaystyle=a_{1}x+a^{3}_{1}\mu
t+b_{1},~{}\eta_{2}=a_{2}x+a^{3}_{2}\mu t+b_{2},$ $\displaystyle a_{1}$
$\displaystyle=1,~{}a_{2}=2,~{}b_{1}=\log(3)/2,~{}b_{2}=b_{1},\mu=1,$

where $\lambda_{1}=1.5$ and $\lambda_{2}=0.25$ are considered.

For the inverse problem, the objective is to identify $\lambda_{1}$ and
$\lambda_{2}$, which are considered unknown, and to predict $u$ in the entire
spatial-temporal domain, given partial noisy measurements on $u$ and $f$.
Specifically, 200 and 100 random noisy measurements on $u$ and $f$,
respectively, are considered available; see also Fig. 7. The noise for the
measurements on both $u$ and $f$ is assumed to be Gaussian with zero mean and
known standard deviation set to 0.1.

Shown in Figs. 7-7 are the errors between the HMC and DEns predicted means,
for $u$ at $t=-1.5$ and 0.4, and the reference solutions bounded by the
predicted total uncertainties. Further, the predicted values for $\lambda_{1}$
and $\lambda_{2}$, which are provided in Table 5, are close to the reference
values for both methods.

| HMC | DEns
---|---|---
$\lambda_{1}$ (mean $\pm$ std) |  1.4095 $\pm$ 0.0372 |  1.5452 $\pm$ 0.0672
$\lambda_{2}$ (mean $\pm$ std) |  0.2642 $\pm$ 0.0057 |  0.3249 $\pm$ 0.0234
Table 5:  Inverse KdV problem: Predictions corresponding to $\lambda_{1}$ and
$\lambda_{2}$ using HMC and DEns; the reference values are $\lambda_{1}=1.5$
and $\lambda_{2}=0.25$.

###  4.4 100-dimensional Darcy problem

Flow through porous media is studied in various disciplines, such as chemical
engineering, e.g., for flows in packed-bed catalytic reactors, oil recovery,
e.g., for the displacement of oil with water, and in geophysics e.g., for the
transport of pollutant in soil. A steady flow through porous media in two
dimensions is considered here, which is described by the Darcy’s law as
follows:

(16) $\displaystyle\nabla\cdot(\lambda(x,y;\xi)\nabla
u(x,y;\xi))=f,~{}x,y\in[0,1],$

where $\lambda$ is the hydraulic conductivity field, $u$ denotes the hydraulic
head, and $f=50$ is considered. The boundary conditions are:

(17a) $\displaystyle u(x=0,y)=1,~{}u(x=1,y)=0,$ (17b)
$\displaystyle\partial_{\bm{n}}u(x,y=0)=\partial_{\bm{n}}u(x,y=1)=0,$

where $\bm{n}$ denotes the unit normal vector of the boundary.

In general, $\lambda(x,y;\xi)$ is determined by the pore structure. In this
study, we employ a stochastic model for $\lambda(x,y;\xi)$ to consider porous
media with different micro-structures. In particular,
$\lambda(x,y;\xi)=\exp(\tilde{\lambda}(x,y;\xi))$, where
$\tilde{\lambda}(x,y;\xi)$ is a truncated Karhunen-Loève expansion of a
Gaussian process with zero mean and the following kernel:

(18a) $\displaystyle
k_{\tilde{\lambda}}(x,y,x^{\prime},y^{\prime})=\exp(-\frac{(x-x^{\prime})^{2}}{2l^{2}}-\frac{(y-y^{\prime})^{2}}{2l^{2}}),$
(18b) $\displaystyle 0\leq x,y,x^{\prime},y^{\prime}\leq 1,~{}l=0.25,$

and we only keep the first 100 leading terms of the expansion.

We consider two different scenarios here, i.e., Case A, in which we employ a
DeepONet mapping $\tilde{\lambda}$ to $u$ pre-trained with 9900 different
paired data ($\tilde{\lambda}$, $u$), and the test data are partial noisy
measurements on $\lambda$ and $u$; and Case B, in which we employ the DEns
method to train a DeepONet using the same training data as in Case A, which
also maps $\tilde{\lambda}$ to $u$. Complete $\lambda$ measurements are used
as test data, i.e., a measurement is available at each required location of
the DeepONet input. Particularly, we refer to the DeepONet trained with DEns
as uncertain DeepONet (U-DeepONet) following [45].

####  4.4.1 Physics-agnostic functional prior

Figure 8:  100-dimensional Darcy problem (Case A): Predictions for
$\log(\lambda)$ and $u$ from PA-GAN-FP+HMC. Training data for $\lambda$ are
shown with circles, whereas training data for $u$ with cross symbols. Total
uncertainty is denoted by $2\bar{\sigma}$.

For Case A, 20 noisy measurements of $\lambda$ and $u$ are assumed to be
available; see also Figs. 8 and 9 for the respective locations. The objective
is to reconstruct the entire $\lambda$ and $u$ given partial noisy
measurements of $\lambda/u$. Specifically, a pre-trained DeepONet combined
with two physics-agnostic models, i.e. a PA-BNN-FP and a PA-GAN-FP, are
considered. In particular, the PA-GAN-FP is pre-trained using the same data
for training the DeepONet to learn the functional prior for $\tilde{\lambda}$,
while we specify a standard normal distribution as the prior for each
parameter in the PA-BNN-FP for the prior of $\tilde{\lambda}$. For posterior
inference, HMC is used for both models. The details for the pre-trained of
DeepONet and PA-GAN-FP can be found in the GitHub repository.

With the pre-trained DeepONet and the physics-agnostic models as the prior in
function space, we present the predictions for $\lambda$ and $u$ in Figs. 8
and 9, respectively. We observe that the errors between the predicted means
and the reference solutions for $\lambda$ and $u$ are mostly bounded by the
predicted uncertainties in both methods. Further, the predicted uncertainties
from PA-GAN-FP are smaller than those from PA-BNN-FP, which is expected since
PA-GAN-FP is equipped with a more informative prior. Interested readers can
refer to [41, 45] for more comparison between the GAN- and BNN-FP.

Figure 9:  100-dimensional Darcy problem (Case A): Predictions for
$\log(\lambda)$ and $u$ from PA-BNN-FP+HMC. Training data for $\lambda$ are
shown with circles, whereas training data for $u$ with cross symbols. Total
uncertainty is denoted by $2\bar{\sigma}$.

####  4.4.2 Uncertain DeepONet

Figure 10:  100-dimensional Darcy problem: Predictions obtained from
U-DeepONet and corresponding to $u$. Total uncertainty is denoted by
$2\bar{\sigma}$.

For Case B, we assume that a complete input $\lambda$ is available, i.e., a
datapoint is available at each required location of $\lambda$. We quantify
model uncertainty in the prediction of $u$ using DEns. In particular, we run
the DeepONet training algorithm five times with different parameter
initializations, and subsequently compute the predicted mean and uncertainties
based on the five predictions of $u$; see Eqs. (7)-(8). The error between the
predicted mean and the reference solution of $u$ for an unseen
$\log({\lambda})$ is bounded by the predicted uncertainties, as shown in Fig.
10.

##  5 Summary

Scientific machine learning (SciML) has emerged recently as an effective and
powerful tool for data fusion, solving ordinary/partial differential equations
(ODEs, PDEs), and learning operator mappings in various scientific and
engineering disciplines. Indicative SciML models are the physics-informed
neural networks (PINNs) for solving forward and inverse ODE/PDE problems and
the deep operator networks (DeepONets) for learning operators that can, for
instance, be used as fast solvers of ODEs/PDEs. In this context, quantifying
predictive uncertainties is crucial for risk-sensitive applications as well as
for efficient and economical design. In this regard, a comprehensive study of
UQ methods for various SciML models, e.g., neural networks (NNs), PINNs,
DeepONets, has been provided in [45].

However, due to the fact that the implementations of current UQ methods are
not straightforward, we have developed in this paper an open-source Python
library (github.com/Crunch-UQ4MI), termed NeuralUQ, that serves as an
efficient and reliable toolbox for UQ in SciML. In this study, we have
presented a detailed tutorial of NeuralUQ, followed by four numerical
examples, including dynamical systems and high-dimensional parametric and
time-dependent PDEs, to demonstrate the applicability, reliability, and
efficiency of NeuralUQ. We have further demonstrated the flexibility and
customizability of NeuralUQ for supporting alternative surrogates, prior and
posterior parameter distributions, as well as the inference methods.
Indicative such modifications have been described and applied in regression
and differential equation problems. More generally, in addition to the UQ
methods presented herein for addressing the problem scenarios of Section 2.1,
NeuralUQ enables the users to address additional scenarios by combining the
built-in components, such as the different surrogates and inference methods.

Despite the plurality of methods implemented in NeuralUQ, appropriately
selecting a UQ method for solving a specific problem is not trivial. In this
regard, we provide a few empirical remarks in the following. Deterministic
methods, such as deep ensemble, are computationally efficient, and they
provide reasonable model uncertainties for both neural ODE/PDE and neural
operators problems. Further, they can handle big data by using data mini-
batches, as in standard neural network (NN) training. For problems with noisy
data, NN parameter regularization can be used for addressing overfitting.
However, this may have a negative effect on predicted uncertainties.
Specifically, a large regularization weight can lead to small model
uncertainty, whereas a small weight may be insufficient for addressing
overfitting issues. As a potential remedy, different regularization weights
can be used and compared by evaluating their corresponding predicted
uncertainties using metrics and validation data. As for Bayesian methods,
Hamiltonian Monte Carlo (HMC) typically achieves both satisfactory
computational accuracy and efficiency for posterior estimation for problems
with relatively small datasets, of the order of a few hundred measurements.
For problems involving big data, the methods Langevin dynamics (LD), mean-
field variational inference (MFVI), and Monte Carlo dropout (MCD), which
support mini-batch training in a straightforward manner, are computationally
more efficient. For solving neural ODE/PDE problems without and with
historical data using Bayesian methods, we recommend the use of a BNN and the
generator of a pre-trained physics-informed generative adversarial network as
surrogates, respectively. Further, as a rule of thumb for neural operator
problems, the physics-agnostic functional prior (PA-FP) can be used to
quantify the uncertainties caused by incomplete data, and the deep ensemble
for quantifying the epistemic uncertainties. Furthermore, the methods HMC and
LD/MFVI are recommended for problems with small and big data, respectively, if
FP is used.

Overall, we envision that NeuralUQ will be used in the future in large-scale
computational science and engineering applications as a toolbox of reliable
and efficient UQ implementations. Clearly, such applications require up-to-
date, well-maintained, and versatile software packages. In this regard,
indicative extensions of our library relate to incorporating more sampling
methods that can handle big data, e.g., stochastic gradient Hamiltonian Monte
Carlo, and variational inference methods with more accurate parameterized
posterior distributions, e.g., normalizing flows. Finally, we will continue
updating the library and we also encourage future users and researchers to
develop additional components for NeuralUQ with the goal of expanding its
current capabilities and advancing the field.

## Acknowledgments

This work is supported by the DOE PhILMs project (No. DE-SC0019453), the MURI-
AFOSR FA9550-20-1-0358 projects, and the DARPA HR00112290029 project.

## References

  * [1] M. Abadi, A. Agarwal, P. Barham, and et. al., TensorFlow: Large-scale machine learning on heterogeneous systems, 2015, https://www.tensorflow.org/.  Software available from tensorflow.org.
  * [2] M. Abdar, F. Pourpanah, S. Hussain, and et. al., A Review of Uncertainty Quantification in Deep Learning: Techniques, Applications and Challenges, arXiv:2011.06225 [cs], (2021), https://arxiv.org/abs/2011.06225.
  * [3] H. Babaee, C. Bastidas, M. DeFilippo, C. Chryssostomidis, and G. Karniadakis, A multifidelity framework and uncertainty quantification for sea surface temperature in the massachusetts and cape cod bays, Earth and Space Science, 7 (2020), p. e2019EA000954.
  * [4] N. Benes, A. Kasman, and K. Young, On decompositions of the KdV 2-soliton, Journal of Nonlinear Science, 16 (2006), pp. 179–200.
  * [5] M. Betancourt, A Conceptual Introduction to Hamiltonian Monte Carlo, arXiv preprint arXiv:1701.02434, (2017), https://arxiv.org/abs/1701.02434.
  * [6] C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra, Weight Uncertainty in Neural Network, in International Conference on Machine Learning, PMLR, 2015, pp. 1613–1622.
  * [7] S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis, Physics-informed neural networks (PINNs) for fluid mechanics: A review, Acta Mechanica Sinica, (2022), pp. 1–12.
  * [8] S. Cai, Z. Wang, F. Fuest, Y. J. Jeon, C. Gray, and G. E. Karniadakis, Flow over an espresso cup: inferring 3-D velocity and pressure fields from tomographic background oriented Schlieren via physics-informed neural networks, Journal of Fluid Mechanics, 915 (2021).
  * [9] G. C. Calafiore, C. Novara, and C. Possieri, A time-varying SIRD model for the COVID-19 contagion in Italy, Annual Reviews in Control, 50 (2020), pp. 361–372.
  * [10] F. Chen, D. Sondak, P. Protopapas, M. Mattheakis, S. Liu, D. Agarwal, and M. Di Giovanni, Neurodiffeq: A python package for solving differential equations with neural networks, Journal of Open Source Software, 5 (2020), p. 1931.
  * [11] Y.-C. Chen, P.-E. Lu, C.-S. Chang, and T.-H. Liu, A time-dependent SIR model for COVID-19 with undetectable infected persons, IEEE Transactions on Network Science and Engineering, 7 (2020), pp. 3279–3294.
  * [12] Y. Chung, I. Char, H. Guo, J. Schneider, and W. Neiswanger, Uncertainty Toolbox: An Open-Source Library for Assessing, Visualizing, and Improving Uncertainty Quantification, arXiv:2109.10254 [cs, stat], (2021), https://arxiv.org/abs/2109.10254.
  * [13] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, and F. Piccialli, Scientific Machine Learning through Physics-Informed Neural Networks: Where we are and What’s next, arXiv preprint arXiv:2201.05624, (2022).
  * [14] J. Darbon, G. P. Langlois, and T. Meng, Overcoming the curse of dimensionality for some Hamilton–Jacobi partial differential equations via neural network architectures, Research in the Mathematical Sciences, 7 (2020), pp. 1–50.
  * [15] R. De Neufville, K. Hodota, J. Sussman, and S. Scholtes, Real options to increase the value of intelligent transportation systems, Transportation Research Record, 2086 (2008), pp. 40–47.
  * [16] J. V. Dillon, I. Langmore, D. Tran, and et. al., Tensorflow distributions, arXiv preprint arXiv:1711.10604, (2017).
  * [17] E. Dong, H. Du, and L. Gardner, An interactive web-based dashboard to track COVID-19 in real time, The Lancet Infectious Diseases, 20 (2020), pp. 533–534.
  * [18] Y. Gal and Z. Ghahramani, Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning, in International Conference on Machine Learning, 2016, pp. 1050–1059.
  * [19] J. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson, Gpytorch: Blackbox matrix-matrix gaussian process inference with GPU acceleration, Advances in Neural Information Processing Systems, 31 (2018).
  * [20] S. Goswami, M. Yin, Y. Yu, and G. E. Karniadakis, A Physics-Informed Variational DeepONet for Predicting the Crack Path in Brittle Materials, arXiv preprint arXiv:2108.06905, (2021), https://arxiv.org/abs/2108.06905.
  * [21] E. Haghighat and R. Juanes, Sciann: A keras/tensorflow wrapper for scientific computations and physics-informed deep learning using artificial neural networks, Computer Methods in Applied Mechanics and Engineering, 373 (2021), p. 113552.
  * [22] J. Han, A. Jentzen, and W. E, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 115 (2018), pp. 8505–8510.
  * [23] O. Hennigh, S. Narasimhan, M. A. Nabian, A. Subramaniam, K. Tangsali, Z. Fang, M. Rietmann, W. Byeon, and S. Choudhry, NVIDIA SimNet™: An AI-Accelerated Multi-Physics Simulation Framework, in International Conference on Computational Science, Springer, 2021, pp. 447–461.
  * [24] M. D. Hoffman, A. Gelman, et al., The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo, J. Mach. Learn. Res., 15 (2014), pp. 1593–1623.
  * [25] G. Huang, Y. Li, G. Pleiss, Z. Liu, J. E. Hopcroft, and K. Q. Weinberger, Snapshot Ensembles: Train 1, get M for free, arXiv:1704.00109 [cs], (2017), https://arxiv.org/abs/1704.00109.
  * [26] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-Informed Machine Learning, Nature Reviews Physics, (2021).
  * [27] K. Kashinath, M. Mustafa, A. Albert, J. Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzadenesheli, R. Wang, A. Chattopadhyay, A. Singh, et al., Physics-informed machine learning: case studies for weather and climate modelling, Philosophical Transactions of the Royal Society A, 379 (2021), p. 20200093.
  * [28] E. Kharazmi, M. Cai, X. Zheng, Z. Zhang, G. Lin, and G. E. Karniadakis, Identifiability and predictability of integer-and fractional-order epidemiological models using physics-informed neural networks, Nature Computational Science, 1 (2021), pp. 744–753.
  * [29] A. Koryagin, R. Khudorozkov, and S. Tsimfer, Pydens: A python framework for solving differential equations with neural networks, arXiv preprint arXiv:1909.11544, (2019).
  * [30] V. Kuleshov, N. Fenner, and S. Ermon, Accurate Uncertainties for Deep Learning Using Calibrated Regression, arXiv preprint arXiv:1807.00263, (2018), https://arxiv.org/abs/1807.00263.
  * [31] B. Lakshminarayanan, A. Pritzel, and C. Blundell, Simple and Scalable Predictive Uncertainty Estimation Using Deep Ensembles, in Advances in Neural Information Processing Systems, 2017, pp. 6402–6413.
  * [32] J. Lao, C. Suter, I. Langmore, C. Chimisov, and et. al., Tfp.mcmc: Modern Markov Chain Monte Carlo Tools Built for Modern Hardware, arXiv:2002.01184 [cs, stat], (2020), https://arxiv.org/abs/2002.01184.
  * [33] D. Levi, L. Gispan, N. Giladi, and E. Fetaya, Evaluating and Calibrating Uncertainty Prediction in Regression Tasks, arXiv preprint arXiv:1905.11659, (2019), https://arxiv.org/abs/1905.11659.
  * [34] G. Lin, C. Moya, and Z. Zhang, Accelerated Replica Exchange Stochastic Gradient Langevin Diffusion Enhanced Bayesian DeepONet for Solving Noisy Parametric PDEs, arXiv preprint arXiv:2111.02484, (2021), https://arxiv.org/abs/2111.02484.
  * [35] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning Nonlinear Operators via DeepONet Based on the Universal Approximation Theorem of Operators, Nature Machine Intelligence, 3 (2021), pp. 218–229.
  * [36] L. Lu, X. Meng, S. Cai, Z. Mao, S. Goswami, Z. Zhang, and G. E. Karniadakis, A Comprehensive and Fair Comparison of Two Neural Operators (with Practical Extensions) Based on FAIR Data, arXiv preprint arXiv:2111.05512, (2021), https://arxiv.org/abs/2111.05512.
  * [37] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, DeepXDE: A deep learning library for solving differential equations, SIAM Review, 63 (2021), pp. 208–228.
  * [38] K. O. Lye, S. Mishra, and D. Ray, Deep learning observables in computational fluid dynamics, Journal of Computational Physics, 410 (2020), p. 109339.
  * [39] X. Meng, H. Babaee, and G. E. Karniadakis, Multi-Fidelity Bayesian Neural Networks: Algorithms and Applications, Journal of Computational Physics, 438 (2021), p. 110361.
  * [40] X. Meng and G. E. Karniadakis, A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems, Journal of Computational Physics, 401 (2020), p. 109020.
  * [41] X. Meng, L. Yang, Z. Mao, J. del Águila Ferrandis, and G. E. Karniadakis, Learning functional priors and posteriors from data and physics, Journal of Computational Physics, 457 (2022), p. 111073.
  * [42] R. M. Neal et al., MCMC using Hamiltonian dynamics, Handbook of Markov Chain Monte Carlo, 2 (2011), p. 2.
  * [43] E. Pickering, G. E. Karniadakis, and T. P. Sapsis, Discovering and forecasting extreme events via active learning in neural operators, arXiv preprint arXiv:2204.02488, (2022).
  * [44] E. Pickering and T. P. Sapsis, Structure and distribution metric for quantifying the quality of uncertainty: Assessing gaussian processes, deep neural nets, and deep neural operators for regression, arXiv preprint arXiv:2203.04515, (2022).
  * [45] A. F. Psaros, X. Meng, Z. Zou, L. Guo, and G. E. Karniadakis, Uncertainty Quantification in Scientific Machine Learning: Methods, Metrics, and Comparisons, arXiv preprint arXiv:2201.07766, (2022).
  * [46] T. Qin, Z. Chen, J. D. Jakeman, and D. Xiu, Deep learning of parameterized equations with applications to uncertainty quantification, International Journal for Uncertainty Quantification, 11 (2021).
  * [47] C. Rackauckas, M. Innes, Y. Ma, J. Bettencourt, L. White, and V. Dixit, Diffeqflux.jl-a julia library for neural differential equations, arXiv preprint arXiv:1902.02376, (2019).
  * [48] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, Universal differential equations for scientific machine learning, arXiv preprint arXiv:2001.04385, (2020), https://arxiv.org/abs/2001.04385.
  * [49] R. Rahaman and A. H. Thiery, Uncertainty Quantification and Deep Ensembles, in Advances in Neural Information Processing Systems, 2021.
  * [50] M. Raissi, P. Perdikaris, and G. 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 (2019), pp. 686–707.
  * [51] M. Raissi, A. Yazdani, and G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science, 367 (2020), pp. 1026–1030.
  * [52] C. Rao, H. Sun, and Y. Liu, Physics informed deep learning for computational elastodynamics without labeled data, arXiv preprint arXiv:2006.08472, (2020).
  * [53] L. Richard de Neufville, S. Scholtes, and T. Wang, Real options by spreadsheet: parking garage case example, Journal of Infrastructure Systems, 12 (2006), pp. 107–111.
  * [54] H. Ritter, A. Botev, and D. Barber, A scalable Laplace approximation for neural networks, in 6th International Conference on Learning Representations, ICLR 2018-Conference Track Proceedings, vol. 6, International Conference on Representation Learning, 2018.
  * [55] G. O. Roberts and J. S. Rosenthal, Optimal scaling of discrete approximations to Langevin diffusions, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60 (1998), pp. 255–268.
  * [56] D. Sen and D. Sen, Use of a modified SIRD model to analyze COVID-19 data, Industrial & Engineering Chemistry Research, 60 (2021), pp. 4251–4260.
  * [57] L. Sun, H. Gao, S. Pan, and J.-X. Wang, Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data, Computer Methods in Applied Mechanics and Engineering, 361 (2020), p. 112732.
  * [58] L. Tierney and J. B. Kadane, Accurate approximations for posterior moments and marginal densities, Journal of the American Statistical Association, 81 (1986), pp. 82–86.
  * [59] R. Tipireddy, P. Stinis, and A. M. Tartakovsky, Time-dependent stochastic basis adaptation for uncertainty quantification, arXiv preprint arXiv:2103.03316, (2021).
  * [60] D. Tran, A. Kucukelbir, A. B. Dieng, M. Rudolph, D. Liang, and D. M. Blei, Edward: A library for probabilistic modeling, inference, and criticism, arXiv preprint arXiv:1610.09787, (2016).
  * [61] N. Trask, A. Huang, and X. Hu, Enforcing exact physics in scientific machine learning: a data-driven exterior calculus on graphs, Journal of Computational Physics, 456 (2022), p. 110969.
  * [62] R. K. Tripathy and I. Bilionis, Deep UQ: Learning deep neural network surrogate models for high dimensional uncertainty quantification, Journal of Computational Physics, 375 (2018), pp. 565–588.
  * [63] X. Wan and G. E. Karniadakis, Multi-element generalized polynomial chaos for arbitrary probability measures, SIAM Journal on Scientific Computing, 28 (2006), pp. 901–928.
  * [64] S. Wang, H. Wang, and P. Perdikaris, Learning the solution operator of parametric partial differential equations with physics-informed DeepONets, Science Advances, 7 (2021), p. eabi8605.
  * [65] M. Welling and Y. W. Teh, Bayesian Learning via Stochastic Gradient Langevin Dynamics, in Proceedings of the 28th International Conference on Machine Learning (ICML-11), 2011, pp. 681–688.
  * [66] J. Willard, X. Jia, S. Xu, M. Steinbach, and V. Kumar, Integrating scientific knowledge with machine learning for engineering and environmental systems, ACM Computing Surveys (CSUR), (2021).
  * [67] A. G. Wilson, The Case for Bayesian Deep Learning, arXiv:2001.10995 [cs, stat], (2020), https://arxiv.org/abs/2001.10995.
  * [68] T. Xifara, C. Sherlock, S. Livingstone, S. Byrne, and M. Girolami, Langevin Diffusions and the Metropolis-adjusted Langevin Algorithm, Statistics & Probability Letters, 91 (2014), pp. 14–19.
  * [69] K. Xu and E. Darve, Adcme: Learning spatially-varying physical fields using deep neural networks, arXiv preprint arXiv:2011.11955, (2020).
  * [70] L. Yang, X. Meng, and G. E. Karniadakis, B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data, Journal of Computational Physics, 425 (2021), p. 109913.
  * [71] L. Yang, D. Zhang, and G. E. Karniadakis, Physics-Informed Generative Adversarial Networks for Stochastic Differential Equations, SIAM Journal on Scientific Computing, 42 (2020), pp. A292–A317.
  * [72] Y. Yang, M. Aziz Bhouri, and P. Perdikaris, Bayesian Differential Programming for Robust Systems Identification under Uncertainty, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476 (2020), p. 20200290.
  * [73] Y. Yang, G. Kissas, and P. Perdikaris, Scalable Uncertainty Quantification for Deep Operator Networks using Randomized Priors, arXiv preprint arXiv:2203.03048, (2022).
  * [74] J. Yao, W. Pan, S. Ghosh, and F. Doshi-Velez, Quality of uncertainty quantification for Bayesian neural network inference, arXiv preprint arXiv:1906.09686, (2019).
  * [75] E. Zelikman, C. Healy, S. Zhou, and A. Avati, CRUDE: Calibrating Regression Uncertainty Distributions Empirically, arXiv:2005.12496 [cs, stat], (2020), https://arxiv.org/abs/2005.12496.
  * [76] D. Zhang, L. Lu, L. Guo, and G. E. Karniadakis, Quantifying Total Uncertainty in Physics-Informed Neural Networks for Solving Forward and Inverse Stochastic Problems, Journal of Computational Physics, 397 (2019), p. 108850.
  * [77] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, and P. Perdikaris, Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data, Journal of Computational Physics, 394 (2019), pp. 56–81.

##  6 Appendix

In this section, we provide a glossary table in Appendix 6.1, and an example
for surrogate customization, i.e., a one-dimensional regression problem, in
Appendix 6.2. We then demonstrate the customization for prior in Samplable and
posterior estimation in Variational, using a one-dimensional differential
equation problem, in Appendix 6.3, and provide the details for data generation
in Appendix 6.4.

###  6.1 Glossary for frequently used terms

In Table 6, we provide a glossary table of frequently used terms in this
paper, and in ML/SciML more generally.

Glossary
---
UQ | uncertainty quantification | SciML | scientific machine learning
DNNs | deep neural networks | BNNs | Bayesian neural networks
GANs | generative adversarial networks | FNN | fully-connected neural network
PINNs | physics-informed neural networks | DeepONets | deep operator networks
FP | functional prior | PI-GANs | physics-informed GANs
PA-FP | physics-agnostic FP | GAN-FP | GAN-based FP
BNN-FP | BNN-based FP |  |
Table 6:  Glossary of frequently used terms in this paper.

###  6.2 Customization of surrogates

We demonstrate the surrogate customization, using a one-dimensional regression
problem. The target function is:

(19) $\displaystyle u=1.5\sin(11x),~{}x\in[-1,1].$

We assume that we have 3 noisy measurements of $u$, which are equidistantly
distributed in $[-0.7,-0.3]$. The noise for the measurements is assumed to be
Gaussian with zero mean and known standard deviation set to 0.05.

We employ three different surrogates for different scenarios: (1) We have
neither historical data nor prior knowledge on the formulation of the target
function. We then employ a BNN with one hidden layer which has 50 neurons and
the hyperbolic tangent activation function as the surrogate. In addition, the
standard normal distribution is utilized as the prior distribution for each
weight and bias in the BNN; (2) the expression describing the target function
is not considered as known, but we have historical data for it. We then train
a GAN to learn the functional prior using historical data; and (3) we have
prior knowledge regarding the target function expressing that it is a sine
function with unknown amplitude and frequency, i.e., $A\sin(\omega x)$, where
$A$ and $\omega$ are parameters to be inferred given measurements. The details
regarding the pre-trained GAN in (2) can be seen in the GitHub repository, and
the posteriors for the parameters in each surrogate are computed using HMC.

Figure 11:  1D regression problem: Predictions for $u$ using different
surrogates. (a) BNN; (b) GAN-FP, and (c) $A\sin(\omega x)$. $2\bar{\sigma}$:
Total uncertainty.

As shown in Fig. 11, the predictions from the GAN-FP and $A\sin({\omega}x)$
surrogates are quite similar for both the predicted means and uncertainties,
which are more accurate than the results from the BNN surrogate in terms of
the means. Also, the errors between the predicted means and the reference
solution are mostly bounded by the predicted uncertainties from GAN-FP and
$A\sin({\omega}x)$, whereas the errors are only partly bounded by the
predicted uncertainties in BNN. Overall, we demonstrate the flexibility and
customizability of NeuralUQ by customizing problem-dependent surrogates to
achieve better predicted accuracy.

###  6.3 Customization of prior and inference

We demonstrate the customization for prior in Samplable and posterior
estimation in Variational, using a one-dimensional differential equation
problem. Specifically, we consider a nonlinear diffusion-reaction system
expressed as follows:

(20) $\displaystyle D\partial_{x}^{2}u-k_{r}u^{3}=f,~{}x\in[-1,1],$

where $u$ represents the concentration of a certain solute, $D$ denotes the
diffusion coefficient, $k_{r}$ is the reaction rate, and $f$ is the source
term.

We first demonstrate the customization of prior using the example of an
inverse problem described by Eq. (20). In this problem, $D=0.01$ is a known
constant, and we have a few noisy measurements on $u$, and $f$. The objective
is to infer $u$ and $f$ for $x\in[-1,1]$, and the reaction rate $k_{r}$ which
is an unknown constant. Specifically, we have noisy measurements of $u$ at 5
random locations, and 17 noisy equidistant measurements of $f$. We then employ
a BNN as well as a Samplable variable as the surrogate for $u$ and $k_{r}$,
respectively. The equation can then be encoded using the automatic
differentiation as in PINNs. The employed BNN has 3 hidden layers with 50
neurons per layers, and the hyperbolic tangent activation function is used as
the activation function. We assume that each weight/bias has the same prior
distribution in this BNN, which is independent Gaussian here, i.e.,
$N(0,1^{2})$. Further, we employ three different prior distributions for
$k_{r}$, i.e., $k_{r}\sim N(0,1^{2})$, $k_{r}\sim HalfN(0.5,0.5^{2})$ and
$\log k_{r}\sim N(0,1^{2})$, where $HalfN$ refers to half-normal distribution.
The posterior samples for the parameters in BNN as well as $k_{r}$ are
obtained using HMC. We only present the predictions for $k_{r}$ in Table 7,
for simplicity. As shown, the predictions for $k_{r}$ are strongly affected by
its prior distribution. Overall, NeuralUQ enables the users to utilize
different prior distributions to achieve better accuracy.

Priors for $k_{r}$ | $k_{r}\sim N(0,1^{2})$ | $k_{r}\sim HalfN(0,1^{2})$ | $\log k_{r}\sim N(0,1^{2})$
---|---|---|---
mean $\pm$ std |  0.5727 $\pm$ 0.4407 |  0.4131 $\pm$ 0.2488 |  0.4538 $\pm$ 0.2550
Table 7:  Nonlinear diffusion-reaction problem: Predictions for $k_{r}$ with
different priors. The reference solution is $k_{r}=0.2$.

Figure 12:  Nonlinear diffusion-reaction problem: Predictions for $u$ and $f$
from VI with different parameterized posteriors. (a) Mean-field Gaussian, (b)
Mean-field Bernoulli.

We proceed to show the customization of posterior estimate in Variational
using the example of a forward problem governed by Eq. (20). In particular,
$D=0.01$, $k_{r}=0.2$, and the boundary condition is $u(-1)=u(1)=0$. In
addition, we assume that we have 10 noisy measurements of $f$ which are
randomly distributed at $x\in[-1,1]$. The noise for the measurements, which is
$N(0,0.01^{2})$, is assumed to be known. The objective is to infer $u$ given
data of $f$. We employ a PI-GAN, which is the same as in [41], to solve this
problem. The same historical data in [41] is utilized here to train the PI-GAN
to obtain the functional priors for $u$ and $f$. Note that HMC is employed in
[41] for the posterior estimation of the parameters in PI-GAN. Here we utilize
VI with two different parameterized posterior distributions to demonstrate the
flexibility of NeuralUQ, i.e., independent Gaussian with unknown means and
standard deviations (mean-field Gaussian), and independent Bernoulli
distributions with fixed $p=0.5$ as the probability of success/failure,
multiplied by unknown variables (mean-field Bernoulli). The predictions for
$u$ and $f$ are displayed in Fig. 12, it is observed that (1) the predicted
means for $u$ and $f$ from the two parameterized posteriors are similar, and
(2) the errors between the predicted means for $u$ and $f$ are mostly bounded
by the predicted uncertainties in both methods.

Finally, note that we employ these specific prior/posterior distributions for
demonstration purposes. Interested users of NeuralUQ can also customize other
problem-dependent prior/posterior distributions by following the examples
presented above.

###  6.4 Data generation

Regarding data generation, in Section 4.1, the nonstiff differential equations
solver ode45 in Matlab is utilized for solving Eq. (9) to generate the
training data as well as the reference solutions; in Section 4.3, the training
data are generated from the exact solution of $u$ from [4]; and in Section
4.4, the finite-element-based Partial Differential Equation Toolbox in Matlab
is utilized for solving Eq. (16) to generate the training data as well as the
reference solutions. Note that all data as well as the employed NNs related to
this study can be accessed in the GitHub repository of NeuralUQ.