# Error-Aware B-PINNs: Improving Uncertainty Quantification in Bayesian
Physics-Informed Neural Networks

Olga Graf
Technical University of Munich
<EMAIL_ADDRESS>Pablo Flores
Pontificia Universidad Católica de Chile
<EMAIL_ADDRESS>
Pavlos Protopapas
Harvard University
<EMAIL_ADDRESS>Karim Pichara
Pontificia Universidad Católica de Chile
<EMAIL_ADDRESS>
###### Abstract

Physics-Informed Neural Networks (PINNs) are gaining popularity as a method
for solving differential equations. While being more feasible in some contexts
than the classical numerical techniques, PINNs still lack credibility. A
remedy for that can be found in Uncertainty Quantification (UQ) which is just
beginning to emerge in the context of PINNs. Assessing how well the trained
PINN complies with imposed differential equation is the key to tackling
uncertainty, yet there is lack of comprehensive methodology for this task. We
propose a framework for UQ in Bayesian PINNs (B-PINNs) that incorporates the
discrepancy between the B-PINN solution and the unknown true solution. We
exploit recent results on error bounds for PINNs on linear dynamical systems
and demonstrate the predictive uncertainty on a class of linear ODEs.

##  1 Introduction

Physics-Informed Neural Networks (PINNs) are deep neural networks (DNNs) that
are able to encode differential equations (DEs) as a component of the neural
network itself. Early results by Lagaris [Lagaris et al., 1998] who first
proposed solving differential equations using neural networks were followed by
a rapid growth of publications on PINNs in the recent years [Raissi et al.,
2019; Hagge et al., 2017; Mattheakis et al., 2019; Mattheakis et al., 2020,
Han et al., 2018; Raissi, 2018]. PINNs have many advantages as compared to
traditional numerical solvers, e.g., they provide closed-form solutions, are
mesh-free. i.e., enable on-demand solution computation after training, can
discover new solutions fast by leveraging transfer learning, and many more.

However, in order to truly flourish and get deployed in safety-critical
applications, deep learning techniques like PINNs need to be _reliable_ ,
i.e., endowed with high-quality uncertainty quantification (UQ) methods. In
the context of PINNs, UQ methods are still scarcely utilized, and there is no
agreed gold standard. Some recent works on the topic include [Psaros et al.,
2022; Zou et al., 2022; Zhang et al., 2018; Graf et al., 2021]. Recently,
physics-informed Gaussian processes (PI-GP) [Raissi et al., 2017] and Bayesian
PINNs (B-PINNs) [Yang et al., 2021] have been proposed as a counterpart to the
classical Bayesian probabilistic machine learning techniques. An alternative
line of research considers the use of generative adversarial networks where
the Bayesian inference is performed in the latent space instead of the
parameter space [Yang et al., 2018; Meng et al., 2021].

UQ in B-PINNs can be seen as more challenging than UQ in traditional numerical
solvers due to multiple sources of uncertainty, ranging from noisy data to
model architecture related issues. There is also a black-box uncertainty
related to how well the network is actually informed about the governing
physics equations. It stems from unsupervised nature of PINNs, but is often
being obscured by considering noisy measurement data. As a rule, all the
forementioned sources of uncertainty are modeled using one "catch-it-all"
term. In this work, we set our goal to disentangle these uncertainties and
show how augmenting separate terms in total uncertainty with useful error
estimates can improve the predictive quality of UQ methods.

Our contributions can be summarized as follows:

  * $\circ$

We empirically and theoretically analyze the aleatoric and the epistemic part
of the total predictive variance in unsupervised B-PINNs. We conclude that
existing B-PINN framework lacks the necessary component to incorporate the
equation-related uncertainty.

  * $\circ$

We introduce the concept of _Error-Aware B-PINNs_ and consider _pseudo-
aleatoric uncertainty_ in the predictive variance.

  * $\circ$

We corroborate our concept in practice on a class of linear dynamical systems,
for which we are able to rigorously model the pseudo-aleatoric uncertainty.

##  2 Background

###  2.1 Problem Formulation

PINNs can solve differential equations expressed, in the most general form, as
follows (cf. [Psaros et al., 2022]):

$\displaystyle\mathcal{F}_{\lambda}[u(x)]=f(x),\quad$ $\displaystyle
x\in\Omega,$ (1) $\displaystyle\mathcal{B}_{\lambda}[u(x)]=b(x),\quad$
$\displaystyle x\in\Gamma,$ (2)

where $x$ is a space-time coordinate in a domain $\Omega$,
$\mathcal{F}_{\lambda}$ is a general differential operator,
$\mathcal{B}_{\lambda}$ is an initial/boundary condition (I/BC) operator
acting on the domain boundary $\Gamma$, $f(x)$ is a source term, $b(x)$ is an
I/BC term, and $\lambda(x)$ are the parameters related to the physics.
Finally, $u(x)$ is the sought solution of the differential equation.

Various problem scenarios for (1)-(2) have been considered in the literature:
$\mathcal{F}_{\lambda}$ and $\mathcal{B}_{\lambda}$ can be either stochastic
or deterministic, known or unknown; $f$, $b$, and $\lambda$ can be either
exactly known or partially known from noisy samples; solution $u$ can be
either completely unknown or partially known from noisy data.

In this paper, we will focus on _data-free deterministic scenario_ with known
and deterministic $\mathcal{F}_{\lambda}$, $\mathcal{B}_{\lambda}$, $f$, $b$,
and $\lambda$. As discussed in the Introduction, this will help us to
disentangle uncertainties arising due to noisy data and/or stochastic
realizations from uncertainties arising due to non-ideal compliance with the
underlying physics.

###  2.2 Physics-Informed Neural Networks

To solve differential equation (1), a DNN, denoted by $u_{\theta}(x)$, is
constructed as an approximation to the true solution $u(x)$. We use $\theta$
to denote all the DNN parameters that will be optimized during training. We
will refer to this network as the _surrogate network_. It is required that
$u_{\theta}(x)$ complies with the physics imposed by the differential
equation. Therefore, we define a _residual network_ ,

$r_{\theta}(x):=\mathcal{F}_{\lambda}[u_{\theta}(x)]-f(x),$ (3)

which shares the same parameters with network $u_{\theta}(x)$ and should
ideally output zero for any input $x\in\Omega$. Network $r_{\theta}(x)$ is an
example of a _Physics-Informed Neural Network_.

There are two main ways of constructing the optimization objective, depending
on how I/BCs are satisfied. So called soft I/BC enforcement incorporates I/BCs
as a penalty term in the loss function. Hard I/BC enforcement encodes I/BCs in
the network design and thus guarantees that they are exactly satisfied. In
this paper, we resort to the latter approach and use a slightly modified
version of Lagaris’ reparameterization technique. E.g., in one-dimensional
case, given an initial condition $u_{0}=u(x_{0})$, we consider a
transformation $\tilde{u}_{\theta}(x)=u_{0}+(1-e^{-(x-x_{0})})u_{\theta}(x)$.
With a slight abuse of notation, hereinafter we refer to $u_{\theta}(x)$ as
the network which has already been reparameterized.

Mean square error (MSE) minimization is widely used in connection with PINNs,
and in case of hard I/BC enforcement it simplifies down to minimization of the
squared residual $r_{\theta}^{2}(x)$. Algorithm 1 summarizes our workflow of
solving a differential equation with unsupervised PINN.

Algorithm 1 Unsupervised PINNs for solving DEs

Input: $\mathcal{F}_{\lambda}$, $\mathcal{B}_{\lambda}$, $f(x)$, $b(x)$,
$\lambda(x)$ 1) Construct $u_{\theta}(x)$ with random initialization of
$\theta$. 2) Transform $u_{\theta}(x)$ to encode
$\mathcal{B}_{\lambda}[u_{\theta}((x)]=b(x)$. 3) Construct
$r_{\theta}(x):=\mathcal{F}_{\lambda}[u_{\theta}(x)]-f(x)$. 4) Train with loss
$r_{\theta}^{2}(x)$ to obtain point estimates $\hat{\theta}$.

Output: $u_{\hat{\theta}}(x)$, $\forall x\in\Omega$

###  2.3 Uncertainty Quantification in Bayesian Neural Networks

Employment of _Bayesian framework_ for UQ in deep learning has proven very
successful [Neal, 1993; MacKay, 1995; Lampinen et al., 2001; Graves, 2011].
Below we briefly review some concepts and methods that will be used throughout
the paper.

Suppose we have a finite set of observations $\mathcal{O}$ on a bounded domain
$\Omega_{\mathcal{O}}\subset\Omega$. Assume that observations are i.i.d. and
drawn from a distribution $p(\mathcal{O}|\mathbf{\theta})$. This distribution
is referred to as _likelihood function_ and is often modeled by a factorizable
Gaussian distribution,

$p(\mathcal{O}|\mathbf{\theta})=\mathcal{N}(\mathcal{O};y_{\theta},\Sigma_{\mathcal{O}}),$
(4)

where the mean is given by a relevant neural network $y_{\theta}$ and
$\Sigma_{\mathcal{O}}=\sigma_{\mathcal{O}}^{2}I$. It is therefore assumed that
observations $\mathcal{O}$ are produced by a data-generating process which
contains deterministic part modeled by $y_{\theta}$ as well as some additive
noise. This noise represents irreducible or _aleatoric uncertainty_ and its
covariance matrix $\Sigma_{\mathcal{O}}$ may be either known or assumed, or
inferred from data.

Standard NN training is performed in a frequentist manner and does not account
for uncertainty. The model parameters $\theta$ are inferred as point estimates
by maximizing the likelihood, which, given that the likelihood is Gaussian, is
equivalent to MSE minimization. Note that in case of unsupervised PINNs, in
order to derive MSE loss as in Algorithm 1, we assume virtual "observations"
which are always equal to zero, whereas the mean is given by a residual
network $r_{\theta}$ and $\sigma_{\mathcal{O}}^{2}$ can be arbitrary.

In contrast to frequentist inference, Bayesian inference places a distribution
not only over the observations $\mathcal{O}$, but over the parameters $\theta$
as well, resulting in Bayesian neural networks (BNNs). The probabilistic
nature of parameters in Bayesian framework gives rise to the _epistemic
uncertainty_.

The _posterior distribution_ over the parameters $\theta$, i.e., distribution
over the "plausible" parameter values based on the observations $\mathcal{O}$,
is given by Bayes’ theorem,

$p(\mathbf{\theta}|\mathcal{O})=\frac{p(\mathcal{O}|\mathbf{\theta})p(\mathbf{\theta})}{p(\mathcal{O})},$
(5)

where $p(\mathbf{\theta})$ is the prior distribution over the parameters, and
$p(\mathcal{O})$ is the evidence. Obtaining the posterior exactly via (5) is
often computationally and analytically intractable. A brief overview of
techniques that can approximate the posterior by obtaining samples from it or
by another distribution is provided in Section 2.3.1.

The goal of UQ is to estimate the _predictive distribution_
$p(\hat{y}|x,\mathcal{O})$. Here, $\hat{y}$ denotes the prediction at $x$, and
is endowed with a conditional distribution on given observations
$\mathcal{O}$. It can be obtained by integrating out the model parameters
$\theta$, or, if intractable, approximated using Monte Carlo estimation,

$p(\hat{y}|x,\mathcal{O})=\int
p(\hat{y}|x,\theta)p(\theta|\mathcal{O})d\theta\approx\frac{1}{N}\sum_{i=1}^{N}p(\hat{y}|x,\theta^{\mathsf{*}}_{i}),$
(6)

where $\theta^{\mathsf{*}}_{i}\sim p(\mathbf{\theta}|\mathcal{O})$ are the
posterior samples.

The mean of predictive distribution can be approximately calculated as

$\mathbb{E}[\hat{y}|x,\mathcal{O}]\approx\frac{1}{N}\sum_{i=1}^{N}y_{\theta^{\mathsf{*}}_{i}}(x):=\tilde{\mu}_{\mathcal{O}}(x).$
(7)

The variance of predictive distribution can also be approximately calculated.
Using the law of total variance and assuming that the distribution
$p(\hat{y}|x,\theta)$ has the same mean and variance as (4), we have

$\displaystyle\mathrm{Var}(\hat{y}|x,\mathcal{O})$
$\displaystyle=\mathbb{E}_{\theta|\mathcal{O}}[\mathrm{Var}(\hat{y}|x,\theta)]+\mathrm{Var}_{\theta|\mathcal{O}}(\mathbb{E}[\hat{y}|x,\theta])$
$\displaystyle=\mathbb{E}_{\theta|\mathcal{O}}[\Sigma_{\mathcal{O}}]+\mathrm{Var}_{\theta|\mathcal{O}}(y_{\theta}(x))$
$\displaystyle\approx\sigma_{\mathcal{O}}^{2}+\frac{1}{N}\sum_{i=1}^{N}(y_{\theta^{\mathsf{*}}_{i}}(x)-\tilde{\mu}_{\mathcal{O}}(x))^{2}:=\widetilde{\sigma}_{\mathcal{O}}^{2}(x).$
(8)

The first and the second summand in (2.3) represent the aleatoric and the
epistemic parts of the total uncertainty, respectively.

####  2.3.1 Methods for Finding the Posterior Distribution

The most popular posterior approximation methods can be divided into two
families: sampling methods and variational methods. There is also an efficient
closed-form solution for the simplified BNN, namely, if we discard parameter
distributions in the hidden layers and keep them only in the last layer. We
will refer to it as Neural Linear Model (NLM).

Sampling Methods. Markov-Chain Monte Carlo methods draw parameter samples from
the posterior by using a Markov chain with $p(\theta|\mathcal{O})$ as its
stationary distribution. Widely used methods are Hamiltonian Monte Carlo and
Langevin Dynamics. The drawback of such methods is the computational cost. We
will not consider them in this paper.

Variational Methods. Variational inference (VI) methods approximate
$p(\theta|\mathcal{O})$ by employing a variational distribution
$q_{w}(\theta)$ that is optimized with respect to $w$ by maximizing the
evidence lower bound. Popular methods are Bayes by Backpropagation ([Blundell
et al., 2015], utilizes factorized Gaussian variational distribution) and
Monte Carlo Dropout ([Srivastava et al., 2014], samples are drawn by randomly
dropping for each sample a fixed percentage of parameters).

Exact Methods. NLM [Snoek et al., 2015] represents compromise between
tractability and expressiveness. One can interpret NLM as a learned feature
basis $\Phi_{w}=[\phi_{w}(x_{1}),...,\phi_{w}(x_{M})]^{T}$ on which Bayesian
linear regression $\Phi_{w}\theta$ is performed. We will denote the residual
network trained with NLM by $r^{\mathsf{NLM}}_{\theta}$ in order to
differentiate from a standard BNN.

Figure 1: Uncertainty estimation with existing B-PINN approach. We solve the
equation $u^{\prime}+3u=f$ with various source terms, using methodolody from
Section 3.1 combined with Variational Inference. We do not aim for the
perfectly trained model on purpose. We see that uncertainty coming from
underfitting is not always captured well, especially in the region where model
was not trained (white region in the figures).

##  3 Uncertainty Quantification in B-PINNs

In the previous section, we have considered all observable information under
one collective term $\mathcal{O}$. In fact, in the context of PINNs, these
observations comprise a variety of sources: noisy and/or limited data,
compliance with the governing physics equations, stochasticity in the physical
model, neural network architecture. We will denote data-related observations
by $\mathcal{D}$, equation-related observations by $\mathcal{P}$, and the
remaining observations collectively by $\mathcal{H}$. As stated earlier, our
primary goal is to study data-free deterministic case with a fixed DNN
architecture, therefore we will focus on modeling the effect of $\mathcal{P}$
on uncertainty. Subsequently, it can be incorporated into the general
framework, resulting in a predictive distribution
$p(\hat{u}|x,\mathcal{D},\mathcal{P},\mathcal{H})$. Here, $\hat{u}$ denotes
the probabilistic prediction for the DE solution.

To the best of our knowledge, equation-related observations for B-PINNs have
not been thoroughly studied in the literature. One reason for that might be
the fact that it is not possible to completely disentangle different sources
of uncertainty:

  * $\circ$

Observations of $f$ can be simultaneously seen as data-related observations
$\mathcal{D}$ and equation-related observations $\mathcal{P}$. Therefore, as
it is not necessary to construct a special likelihood function for
$\mathcal{P}$, this might obscure the importance of taking $\mathcal{P}$
explicitly into account at the later inference stages.

  * $\circ$

The total variance in (2.3) is affected by aleatoric uncertainty not only via
the corresponding term, but also via the epistemic term, which is expected to
increase with more noise. This can obscure the importance of properly modeling
the aleatoric term $\sigma_{\mathcal{P}}^{2}$ in
$\mathrm{Var}(\hat{u}|x,\mathcal{P})$. In fact, this term is absent in the
existing B-PINN methodology, which might result in less reliable uncertainty
estimates.

###  3.1 Existing Approach

In this subsection, we will apply and analyze the original approach by Yang et
al. [Yang et al., 2021] for B-PINNs in data-free case. We will empirically
show that it does not always lead to useful uncertainty estimates, especially
in the region where model was not trained.

Since we assume that $f$ is exactly known and there is no noise, we can
consider a virtual dataset $\mathcal{D}=\\{(x_{j},0)\\}_{j=1}^{M}$, whereas
the corresponding likelihood is given by

$p(\mathcal{D}|\mathbf{\theta})=\prod_{j=1}^{M}p(0|x_{j},\mathbf{\theta})=\prod_{j=1}^{M}\mathcal{N}(0;r_{\theta}(x_{j}),\sigma_{\mathcal{D}}^{2}).$
(9)

If $\sigma_{\mathcal{D}}^{2}$ is not known, it can be estimated, e.g., via
hierarchical modeling. The prior distribution $p(\theta)$ is assumed to be
Gaussian with zero mean and diagonal covariance matrix. Corresponding
variances, if not known, can be inferred similarly to
$\sigma_{\mathcal{D}}^{2}$. Then, after we obtain samples from the posterior
using methodology from Section 2.3.1, we estimate the predictive distribution

$p(\hat{u}|x,\mathcal{D})=\mathbb{E}_{\theta|\mathcal{D}}[\hat{u}|x,\theta].$
(10)

Since we are in data-free scenario, a natural choice for the distribution on
the right-hand side of (10) is

$p(\hat{u}|x,\theta)=\delta(\hat{u}-u_{\theta}(x)).$ (11)

Substituting (11) into (2.3) yields

$\displaystyle\mathrm{Var}(\hat{u}|x,\mathcal{D})\approx\widetilde{\sigma}_{\mathcal{D}}^{2}(x)=\frac{1}{N}\sum_{i=1}^{N}(u_{\theta^{\mathsf{*}}_{i}}(x)-\tilde{\mu}_{\mathcal{D}}(x))^{2},$
(12)

where only the epistemic term is present.

Although the predictive variance in (12) is "informed" about the underlying DE
implicitly via $\mathcal{D}$, it might not be large enough to cover the true
solution outside of the training region or when B-PINN is for some reason not
trained well, as shown in Figure 1.

Another drawback of the current approach is its incompatibility with
$r^{\mathsf{NLM}}_{\theta}$. By imposing
$\mathcal{F}_{\lambda}[\Phi_{w}\theta]$, we lose the ability to derive a
closed-form solution.

###  3.2 Error-Aware B-PINNs

In classical function approximation problem with BNN, there is no guarantee
that uncertainty covers the error $y-y_{\theta}$ in the region where no data
is available, unless there are no specific assumptions on the function class.
Unlike BNNs, B-PINNs are to a certain extent aware of discrepancy between the
surrogate network and the true solution. Of course, this depends on the choice
of the loss function as well as on the link between the chosen loss function
and the error that the surrogate network makes.

In this subsection, we introduce a _pseudo-aleatoric uncertainty_ that
quantifies the discrepancy between $u$ and $u_{\theta}$. We still call it
aleatoric because it acts as a counterpart of aleatoric uncertainty in the
classical Bayesian framework. On the other hand, it is not coming from data,
but from PINN itself.

Our goal is to construct a distribution $p(\hat{u}|x,\theta,\mathcal{P})$ with
the mean given by $u_{\theta}(x)$ and the pseudo-aleatoric variance given by
some operator $\mathcal{E}$ which acts on pre-trained PINN,

$\sigma_{\mathcal{P}}(x)=\mathcal{E}[r_{\mathsf{MSE}}(x),\lambda(x)],$ (13)

where $r_{\mathsf{MSE}}$ denotes PINN trained with MSE. Clearly, this is not a
straightforward task which raises questions how to construct $\mathcal{E}$ and
how to find an appropriate functional form of the PDF. At this point, we are
able to answer the first question with provable guarantees in the special case
of linear dynamical systems (see Section 3.2.1). As for functional form of
$p(\hat{u}|x,\theta,\mathcal{P})$, we consider normal distribution as a proof
of concept. Although it does not represent the probability of the true
solution, it makes sure that the uncertainty is inflated in the region where
the error is high and that the true solution is covered by uncertainty.

The predictive variance in the error-aware case is given by

$\displaystyle\mathrm{Var}(\hat{u}|x,\mathcal{D},\mathcal{P})$
$\displaystyle\approx\widetilde{\sigma}_{\mathcal{P}\mathcal{D}}^{2}(x)=\sigma_{\mathcal{P}}^{2}(x)+\frac{1}{N}\sum_{i=1}^{N}(u_{\theta^{\mathsf{*}}_{i}}(x)-\tilde{\mu}_{\mathcal{D}}(x))^{2}.$
(14)

In case of DEs that do not fall into category of Section 3.2.1, there is no
rigorous theory on how to construct (13). Various heuristic approaches can be
considered, whereas the simplest one is equating
$\mathcal{E}[r_{\mathsf{MSE}}(x)]=r_{\mathsf{MSE}}(x)$ or
$\mathcal{E}[r_{\mathsf{MSE}}(x_{n})]=\sum_{i=0}^{n}r_{\mathsf{MSE}}(x_{i})$.
We test the latter in Section 3.2.2 on an example of a non-linear PDE.

Now we show how to combine our approach with $r^{\mathsf{NLM}}_{\theta}$ such
that the closed-form solution still exists. Namely, our distribution with
pseudo-aleatoric variance can be used as a likelihood function for simulated
data
$\mathcal{D}_{\mathcal{P}}=\\{(x_{j},u_{\mathsf{MSE}}(x_{j}))\\}_{j=1}^{M}$,
where $u_{\mathsf{MSE}}$ denotes surrogate network trained with MSE, resulting
in

$\displaystyle p(\mathcal{D}_{\mathcal{P}}|\mathbf{\theta})$
$\displaystyle=\prod_{j=1}^{M}\mathcal{N}(u_{\mathsf{MSE}}(x_{j});u_{\theta}(x_{j}),\sigma_{\mathcal{P}}^{2}(x_{j})).$
(15)

If used in conjunction with this likelihood, $r^{\mathsf{NLM}}_{\theta}$
provides tractable inference. Note that $\sigma_{\mathcal{P}}^{2}(x_{j})$ is
now heteroscedastic. As long as Bayesian linear regression model with
heteroscedastic variance is less commonly used, we give here an explicit form
of the predictive distribution,

$\displaystyle p(\hat{u}_{\mathsf{NLM}}|x,\mathcal{D}_{\mathcal{P}})$
$\displaystyle=\mathcal{N}(\hat{u}_{\mathsf{NLM}};\mu_{\mathcal{D}_{\mathcal{P}}}(x),\sigma^{2}_{\mathcal{D}_{\mathcal{P}}}(x))$
(16)
$\displaystyle=\mathcal{N}(\hat{u}_{\mathsf{NLM}};\Phi_{w}\mu_{\mathsf{p}},\sigma^{2}_{\mathcal{P}}(x)+\Phi_{w}\Sigma_{\mathsf{p}}\Phi_{w}^{T}),$

where posterior covariance matrix and mean are given by

$\displaystyle\Sigma_{\mathsf{p}}$
$\displaystyle=(\Phi_{w}^{T}\Sigma_{\mathcal{P}}\Phi_{w}+\sigma^{-2}I)^{-1},$
(17) $\displaystyle\mu_{\mathsf{p}}$
$\displaystyle=\Sigma_{\mathsf{p}}(\Sigma_{\mathcal{P}}^{-1}\Phi_{w})^{T}u_{\mathsf{MSE}}(x_{\mathcal{D}}),$
(18)

$x_{\mathcal{D}}=[x_{1},...,x_{M}]^{T}$,
$\Sigma_{\mathcal{P}}=\mathrm{diag}([\sigma_{\mathcal{P}}^{2}(x_{1}),...,\sigma_{\mathcal{P}}^{2}(x_{M})])$,
and $\sigma^{2}$ denotes prior variance.

Moreover, this setup for $r^{\mathsf{NLM}}_{\theta}$ allows for prior
optimization. We require the predictive uncertainty to cover the error bounds
given by $\sigma_{\mathcal{P}}(x)$, i.e.,

$\displaystyle|\mu_{\mathcal{D}_{\mathcal{P}}}(x)-u_{\mathsf{MSE}}(x)|\leqslant
3\sigma_{\mathcal{D}_{\mathcal{P}}}(x)-\sigma_{\mathcal{P}}(x).$ (19)

We only consider priors that satisfy the above inequality. To find $\sigma$
that makes the posterior predictive follow the error bounds as closely as
possible, we solve the following minimization problem,

$\displaystyle\underset{\sigma}{\min}\|\mu_{\mathcal{D}_{\mathcal{P}}}-u_{\mathsf{MSE}}\|+\|\sigma_{\mathcal{D}_{\mathcal{P}}}-\sigma_{\mathcal{P}}\|.$
(20)

This simple approach can be regarded as error visualization technique. It can
be implemented on top of deterministically trained network at almost no
computational cost, yet provides quick and reasonable inference for the
goodness of fit. We show the results for $r^{\mathsf{NLM}}_{\theta}$ with
optimal prior in Figure 2.

Figure 2: Uncertainty estimation with Error-Aware B-PINNs. We solve the
equation $u^{\prime}+3u=f$ with various source terms, using methodolody from
Section 3.2 combined with Variational Inference (VI-B-PINN) and Neural Linear
Model (NLM-PINN). Our approach shows reasonable uncertainty estimates even in
the test region where the model was not trained (white region in the figures)
and in case it was trained only for 10 epochs. Figure 3: Uncertainty
estimation for Burgers’ equation. We compare the existing B-PINN method (upper
row) with our method (lower row). Uncertainty coming from underfitting is also
not captured well due to our heuristic approach, nevertheless, provides more
reasonable estimates than in the first case.

Finally, we summarize our method for Error-Aware B-PINNs in Algorithm 2.

Algorithm 2 Error-Aware B-PINNs for solving DEs

Input: $\mathcal{F}_{\lambda}$, $\mathcal{B}_{\lambda}$, $f(x)$, $b(x)$,
$\lambda(x)$ 1) Construct and train $r_{\mathsf{MSE}}(x)$ as in Algorithm 1.
2) Use $r_{\mathsf{MSE}}(x)$ to construct pseudo-aleatoric uncertainty term
$\sigma_{\mathcal{P}}(x)$. 3) Construct Bayesian model $u_{\theta}(x)$ with
$p(\theta)$. 4) Transform $u_{\theta}(x)$ to encode
$\mathcal{B}_{\lambda}[u_{\theta}((x)]=b(x)$. 5) Construct B-PINN
$r_{\theta}(x):=\mathcal{F}_{\lambda}[u_{\theta}(x)]-f(x)$. 6) Obtain the
posterior distribution of $\theta$ via one of the techniques from Section
2.3.1. Likelihood (9) or (15) can be used for this purpose. 7) Estimate mean
and variance of $p(\hat{u}|x,\mathcal{D},\mathcal{P})$ via (3.2).

Output: Estimated $p(\hat{u}|x,\mathcal{D},\mathcal{P})$, $\forall x\in\Omega$

####  3.2.1 Motivating Example: Error-Aware B-PINNs for Linear Dynamical
Systems

One way to encode discrepancy between the surrogate network and the true
solution into the distribution of $\hat{u}$ is to encorporate error bounds
$|u(x)-u_{\theta}(x)|$ in its variance. We emphasize that while this does
_not_ give the probability of seeing the true solution, it is a good starting
point to show the importance of the aleatoric term in the total variance.

Recently, there appeared several works that study the failure modes and
absolute error of PINN solutions [Liu et al., 2021; Graf et al, 2021]. In this
subsection, we employ the approach of Liu et al. [Liu et al., 2021], which
considers unsupervised scenario with hard I/BCs and derives error bounds which
only depend on $x$ and $r_{\theta}$. It makes no assumption on the network
architecture or whether the network is suffciently trained.

E.g., for the first order linear ODE,

$\displaystyle u^{\prime}(x)+\lambda u(x)=f(x),\quad\lambda>0,\quad
u(x_{0})\neq 0,$ (21)

the error bound is given by

$\displaystyle|u(x)-u_{\theta}(x)|\leqslant e^{-\lambda
x}\int_{x_{0}}^{x}e^{\lambda\xi}|r_{\theta}(\xi)|d\xi,$ (22)

provided that $|r_{\theta}|$ is bounded on $\Omega$. We note that in practice
we approximate this integral by summation, as we cannot evaluate the residual
at every point in $\mathbb{R}^{n}$.

Inspired by this bound, we take

$\sigma_{\mathcal{P}}(x)=\mathcal{E}[r_{\mathsf{MSE}}(x),\lambda]=e^{-\lambda
x}\int_{x_{0}}^{x}e^{\lambda\xi}|r_{\mathsf{MSE}}(\xi)|d\xi.$ (23)

Assuming normal distribution, we obtain the results for several source terms
and two inference strategies, Variational Inference (VI-B-PINN) and Neural
Linear Model (NLM-PINN), as shown in Figure 2.

Existing methodolody allows to consider also higher order linear ODE with
constant and non-constant $\lambda(x)$, and systems of first order linear ODEs
with constant coefficients. We present the results for the second order ODEs
in Appendix.

####  3.2.2 Outlook: Error-Aware B-PINNs for PDEs

We consider Burgers’ equation (non-linear PDE),

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=\nu\frac{\partial^{2}u}{\partial x^{2}},$ (24)

as a proof of concept that our methodology might prove useful also in more
advanced cases. The results are presented in Figure 3.

##  4 Conclusion and Future Work

Uncertainty quantification is a challenging task for deep learning models,
especially if it is critical that the predictive uncertainty reflects the
phenomena under consideration and is not giving misleading results.

In this work, we tried to improve the UQ in Bayesian Physics-Informed Neural
Networks by enforcing the information about network’s compliance with the
underlying differential equation. For that purpose, we introduced Error-Aware
B-PINNs with the extra term dubbed as pseudo-aleatoric variance. We showed how
to obtain this variance using the residual-based error bounds in case of
linear dynamical systems.

As a future extension, several non-linear ODEs can be considered by exploiting
solution via perturbation series and constructing similar error bounds and
error bound based variances. A more ambitious approach would be to establish
distributions for the solution errors as opposed to the merely qualitative
assessment of the presented approach.

## References

  * Blundell et al., 2015 Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D. (2015).  Weight uncertainty in neural networks.  In International Conference on Machine Learning, pages 1613–1622.
  * Graf et al., 2021 Graf, O., Flores, P., Protopapas, P., and Pichara, K. (2021).  Uncertainty quantification in neural differential equations.  arXiv preprint arXiv:2111.04207.
  * Graves, 2011 Graves, A. (2011).  Practical variational inference for neural networks.  In Advances in Neural Information Processing Systems, pages 2348–2356.
  * Hagge et al., 2017 Hagge, T., Stinis, P., Yeung, E., and Tartakovsky, A. M. (2017).  Solving differential equations with unknown constitutive relations as recurrent neural networks.  arXiv preprint arXiv:1710.02242.
  * Han et al., 2018 Han, J., Jentzen, A., and E, W. (2018).  Solving high-dimensional partial differential equations using deep learning.  Proceedings of the National Academy of Sciences, 115(34):8505–8510.
  * Lagaris et al., 1998 Lagaris, I., Likas, A., and Fotiadis, D. (1998).  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000.
  * Lampinen and Vehtari, 2001 Lampinen, J. and Vehtari, A. (2001).  Bayesian approach for neural networks – review and case studies.  Neural Networks, 14:257–274.
  * Liu et al., 2021 Liu, S., Huang, X., and Protopapas, P. (2021).  Evaluating error bound for physics-informed neural networks on linear dynamical systems.  arXiv preprint arXiv:2207.01114.
  * MacKay, 1995 MacKay, D. J. (1995).  Bayesian methods for neural networks: Theory and applications.  Neural Networks Summer School, University of Cambridge.
  * Mattheakis et al., 2019 Mattheakis, M., Protopapas, P., Sondak, D., Giovanni, M. D., and Kaxiras, E. (2019).  Physical symmetries embedded in neural networks.  arXiv preprint arXiv:1904.08991.
  * Mattheakis et al., 2020 Mattheakis, M., Sondak, D., Dogra, A. S., and Protopapas, P. (2020).  Hamiltonian neural networks for solving differential equations.  arXiv preprint arXiv:2001.11107.
  * Meng et al., 2021 Meng, X., Yang, L., Mao, Z., del Aguila Ferrandis, J., and Karniadakis, G. E. (2021).  Learning functional priors and posteriors from data and physics.  Journal of Computational Physics, 457:111073.
  * Neal, 1993 Neal, R. M. (1993).  Probabilistic inference using Markov chain Monte Carlo methods.  Department of Computer Science, University of Toronto, Toronto, Ontario, Canada.
  * Psaros et al., 2022 Psaros, A. F., Meng, X., Zou, Z., Guo, L., and Karniadakis, G. E. (2022).  Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons.  arXiv preprint arXiv:2201.07766.
  * Raissi, 2018 Raissi, M. (2018).  Forward-backward stochastic neural networks: Deep learning of high-dimensional partial differential equations.  arXiv preprint arXiv:1804.07010.
  * Raissi et al., 2017 Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2017).  Machine learning of linear differential equations using Gaussian processes.  Journal of Computational Physics, 348:683–693.
  * Raissi et al., 2019 Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2019).  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.
  * Snoek et al., 2015 Snoek, J., Rippel, O., Swersky, K., Kiros, R., Satish, N., Sundaram, N., Patwary, M. M. A., Prabhat, and Adams, R. P. (2015).  Scalable bayesian optimization using deep neural networks.  In International Conference on Machine Learning.
  * Srivastava et al., 2014 Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014).  Dropout: A simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research, 15(56):1929–1958.
  * Yang et al., 2021 Yang, L., Meng, X., and Karniadakis, G. E. (2021).  B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data.  Journal of Computational Physics, 425:109913.
  * Yang and Perdikaris, 2018 Yang, Y. and Perdikaris, P. (2018).  Adversarial uncertainty quantification in physics-informed neural networks.  Journal of Computational Physics, 394:136–152.
  * Zhang et al., 2018 Zhang, D., Lu, L., Guo, L., and Karniadakis, G. E. (2018).  Quantifying total uncertainty in physics-informed neural networks for solving forward and inverse stochastic problems.  Journal of Computational Physics, 397:108850.
  * Zou et al., 2022 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.