remarkRemark hypothesisHypothesis claimClaim L-HYDRAZ. Zou and G. E.
Karniadakis

# L-HYDRA: Multi-Head Physics-Informed Neural Networks

Zongren Zou Division of Applied Mathematics, Brown University, Providence, RI
02912, USA<EMAIL_ADDRESS>George Em Karniadakis Corresponding
author. Division of Applied Mathematics, Brown University, Providence, RI
02912, USA<EMAIL_ADDRESS>
###### Abstract

We introduce multi-head neural networks (MH-NNs) to physics-informed machine
learning, which is a type of neural networks (NNs) with all nonlinear hidden
layers as the body and multiple linear output layers as multi-head. Hence, we
construct multi-head physics-informed neural networks (MH-PINNs) as a potent
tool for multi-task learning (MTL), generative modeling, and few-shot learning
for diverse problems in scientific machine learning (SciML). MH-PINNs connect
multiple functions/tasks via a shared body as the basis functions as well as a
shared distribution for the head. The former is accomplished by solving
multiple tasks with MH-PINNs with each head independently corresponding to
each task, while the latter by employing normalizing flows (NFs) for density
estimate and generative modeling. To this end, our method is a two-stage
method, and both stages can be tackled with standard deep learning tools of
NNs, enabling easy implementation in practice. MH-PINNs can be used for
various purposes, such as approximating stochastic processes, solving multiple
tasks synergistically, providing informative prior knowledge for downstream
few-shot learning tasks such as meta-learning and transfer learning, learning
representative basis functions, and uncertainty quantification. We demonstrate
the effectiveness of MH-PINNs in five benchmarks, investigating also the
possibility of synergistic learning in regression analysis. We name the open-
source code “Lernaean Hydra” (L-HYDRA), since this mythical creature possessed
many heads for performing important multiple tasks, as in the proposed method.

###### keywords:

PINNs, meta-learning, multi-tasking, transfer learning, generative models,
normalizing flows, stochastic problems

34F05, 62M45, 65L99, 65M99, 65N99

##  1 Introduction

Learning across tasks has drawn great attention recently in deep learning and
is an emerging theme in scientific machine learning (SciML), due to the fact
that several classes of scientific problems are similar and/or related
intrinsically by their common physics. Intuitively, if tasks are similar,
e.g., in the context of approximating stochastic processes [44], learning
solution operators of ordinary/partial differential equations (ODEs/PDEs)
[28], and solving parametric PDEs [42, 19, 4], it may be beneficial to relate
them in the modeling, algorithm design, and/or solving procedure. In this
regard, machine learning solvers, developed rapidly in the past few years, are
considerably more flexible and of higher potential compared to traditional
numerical solvers. Significant progress has been witnessed in the general
area, including meta-learning for solving ODEs/PDEs [30, 27, 34, 6], transfer
learning for physics-informed neural networks (PINNs) [3, 7], transfer
learning for domain shift in solving PDEs [14], multi-task learning for PINNs
[40], and generative methods for solving stochastic differential equations
(SDEs) [44, 46, 15]. More recently, operator learning [28, 24] in which direct
operator mapping is learned and subsequently used for other tasks in one-shot
format has attracted a lot of attention.

Multi-head neural networks (MH-NNs) fit perfectly different scenarios of
learning across tasks. They were originally proposed as members of hard-
parameter sharing neural networks (NNs) for deep multi-task learning (MTL)
[5], in which multiple tasks, denoted as $\mathcal{T}_{k},k=1,...,M$, where
$M$ is the number of total tasks, are solved simultaneously. The general goals
of using MH-NNs in MTL are diverse: achieving better performance for all
tasks, learning good and useful representations for downstream tasks, and/or
boosting the learning of main tasks with the help of auxiliary tasks.
Moreover, although originally designed for solving multiple tasks, MH-NNs in
recent years have also been extensively used for meta-learning. For example,
in [41], the connection between MTL and meta-learning was analyzed, and meta-
learning algorithms for MH-NN were discussed; in [25], it was shown that MH-
NNs, trained in MTL fashion also perform task-specific adaptation in meta-
learning; [37] argued that the effectiveness of model-agnostic meta-learning
[10], a well-known meta-learning algorithm, may be due to successfully learned
good representations rather than learned adaptation, and MH-NNs were used to
study the detailed contributions of NNs in fast task adaptations. Overall, it
is commonly acknowledged in the literature that when used to solve previous
tasks, MH-NNs are capable of distilling useful shared information and storing
it in their bodies and heads.

In this paper, we develop MH-NNs for physics-informed machine learning [17],
propose multi-head physics-informed neural networks (MH-PINNs), and further
investigate their applicability and capabilities to MTL, generative modeling,
and meta-learning. A MH-PINN, as shown in Fig. 1, is built upon a conventional
MH-NN and consists of two main parts, the body and multiple heads, and each
head connects to a specific ODE/PDE task. Many architecture splitting
strategies for MH-NNs are adopted in different applications scenarios; e.g.,
for some computer vision problems, a NN is split such that the body consists
of convolutional layers and is followed by fully-connected layers as heads. In
this paper, however, we choose the simplest one, i.e., the body consists of
all nonlinear layers and the head is the last linear layer, for the following
two reasons: (1) the dimensionality of the head is reduced, which enables fast
density estimation (see next section); and (2) the body spontaneously provides
a set of basis functions.

Figure 1: Schematic view of the structure of multi-head physics-informed
neural networks (MH-PINNs) with $M$ different heads, which are built upon
conventional multi-head neural networks. The shared layers are often referred
to as body and the task-specific layer as head. Generally, $u_{k},k=1,...,M$
represent $M$ solutions to $M$ different ODEs/PDEs, formulated in Eq. (1),
which may differ in source terms $f_{k}$, boundary/initial condition terms
$b_{k}$, or differential operator $\mathcal{F}_{k}$.

The novelty and major contributions of this work are as follows:

  1. 1.

We propose a new physics-informed generative method using MH-PINNs for
learning stochastic processes from data and physics.

  2. 2.

We propose a new method for physics-informed few-shot regression problems with
uncertainty quantification using MH-PINNs.

  3. 3.

We study and demonstrate the effectiveness of MTL and synergistic learning
with MH-NNs in regression problems.

The paper is organized as follows. In Sec. 2, we present the problem
formulation, details of MH-PINNs, and the general methodology, including how
to use MH-PINNs for MTL, generative modeling, downstream few-shot physics-
informed learning with uncertainty quantification (UQ). In Sec. 3, we discuss
existing research closely related to our work and compare them conceptually.
In Sec. 4, we test MH-PINNs with five benchmarks, each of which corresponds to
one or more learning purposes, e.g., MTL and generative modeling. In Sec. 5,
we investigate MTL and synergistic learning with the function approximation
example. We conclude and summarize in Sec. 6. The details of our experiments,
such as NN architectures and training strategies, can be found in Appendix A
and B, as well as in the L-HYDRA open-source codes on GitHub, which will be
released once the paper is accepted.

##  2 Methodology

We assume that we have a family of tasks, $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$,
each of which is associated with data $\mathcal{D}_{k},k=1,...,M$. The primary
focus of this paper is on scientific computing and ODEs/PDEs, and therefore we
further assume $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ are physics-informed
regression problems [17].

Consider a PDE of the following form:

(1a) $\displaystyle\mathcal{F}_{k}[u_{k}(x)]$
$\displaystyle=f_{k}(x),x\in\Omega_{k},$ (1b)
$\displaystyle\mathcal{B}_{k}[u_{k}(x)]$
$\displaystyle=b_{k}(x),x\in\partial\Omega_{k},$

where $k$ denotes the index of the task and $k=1,...,M$, $x$ is the general
spatial-temporal coordinate of $D_{x}$ dimensions, $\Omega_{k}$ are bounded
domains, $f_{k}$ and $u_{k}$ are the $D_{u}$-dimensional source terms and
solutions to the PDE, respectively, $\mathcal{F}_{k}$ are general differential
operators, $\mathcal{B}_{k}$ are general boundary/initial condition operators,
and $b_{k}$ are boundary/initial condition terms. For simplicity, throughout
this paper, the domain and the boundary/initial operator, denoted as $\Omega$
and $\mathcal{B}$, are assumed to be the same for all tasks, and the solutions
$u_{k}$ to be task-specific. The task $\mathcal{T}_{k}$ is described as
approximating $u_{k}$, and/or $f_{k}$, and/or $\mathcal{F}_{k}$, and/or
$b_{k}$, from data $\mathcal{D}_{k}$ and Eq. (1).

Traditional numerical solvers often tackle $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$
independently, without leveraging or transferring knowledge across tasks. The
PINN method [38] was designed to solve ODEs/PDEs independently using NNs,
which, however, yields $M$ uncorrelated results. In this paper instead we
treat $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ as a whole and connect them with MH-
PINNs, the architecture of which, shown in Fig. 1, enforces basis-functions-
sharing predictions on the solutions $u_{k}$. In addition to the informative
representation/body, we further relate $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ by
assuming that their corresponding heads in MH-PINNs, denoted as
$\\{H_{k}\\}_{k=1}^{M}$, are samples of a random variable with unknown
probability density function (PDF), denoted as $H$ and $p(H)$, respectively.
The shared body and a generative model of $H$ immediately form a generative
model of the solution $u$, and generators of the source term $f$ and the
boundary/initial term $b$ as well by substituting $u$ into Eq. (1) and
automatic differentiation [1], from which a generative method for
approximating stochastic processes is seamlessly developed.

Generators of $u$, $f$ and $b$, as discussed in [30], are able to provide an
informative prior distribution in physics-informed Bayesian inference [43, 26]
as well as in UQ for SciML [47, 36], where the informative prior compensates
for the insufficiency of observational data to address the physics-informed
learning problems with even a few noisy measurements. In this paper, we
generalize such problem to deterministic cases as well, where the data is
noiseless and methods and results are deterministic, and refer to it as few-
shot physics-informed learning. The general idea is to apply prior knowledge
learned from connecting $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ with MH-PINNs to new
tasks, denoted as $\tilde{\mathcal{T}}$, associated with insufficient data
$\tilde{\mathcal{D}}$, for accurate and trustworthy predictions. The schematic
view of the learning framework is illustrated in Fig. 2, and the details are
explained next.

Figure 2: Schematic view of the learning framework and the proposed method.
Three general types of learning are addressed: physics-informed learning,
generative modeling, and few-shot learning. The physics-informed learning is
performed with MH-PINNs; the generative modeling is done afterwards by density
estimate over the head via normalizing flows (NFs); in the end the few-shot
physics-informed learning is accomplished with prior knowledge obtained from
previous two via either fine-tuning with the learned regularization or
Bayesian inference with the learned prior distribution. The body represents
the set of basis functions learned from solving
$\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ with MH-PINNs, and the density of the head,
estimated from its samples using NFs, acts as the regularization, the prior
distribution, or the generator together with the body, depending on the usage
of MH-PINNs in applications.

###  2.1 Multi-head physics-informed neural networks (MH-PINNs)

Hard parameter sharing is the most commonly used approach when MTL with NNs
are considered, and MH-NNs, as its simplest instance, are frequently adopted
[5, 39]. A MH-PINN, as described earlier, is composed of a body and multiple
heads. We denote by $\Phi$ the body and by $H_{k}$ the head for
$\mathcal{T}_{k}$. Notice that here
$\Phi:\mathbb{R}^{D_{x}}\rightarrow\mathbb{R}^{L}$ is a function parameterized
by a neural network with parameter $\theta$, and $H_{k}\in\mathbb{R}^{L+1}$ is
a vector, where $L$ is the number of neurons on the last layer of the body.
Let us define $H_{k}=[h_{k}^{0},h_{k}^{1},...,h_{k}^{L}]^{T}$,
$\Phi(x)=[\phi^{1}(x),...,\phi^{L}(x)]^{T}$, where
$\phi:\mathbb{R}^{D_{x}}\rightarrow\mathbb{R}$, and then the surrogate for the
solution in $\mathcal{T}_{k}$ can be rewritten as
$\hat{u}_{k}(x)=h_{k}^{0}+\sum_{l=1}^{L}h^{l}_{k}\phi^{l}(x),\forall
x\in\Omega$. The approximated source terms and boundary/initial terms are
derived from Eq. (1) accordingly. In the MTL framework , given data
$\\{\mathcal{D}_{k}\\}_{k=1}^{M}$ and physics Eq. (1), the loss function
$\mathcal{L}$ is formulated as follows:

(2)
$\mathcal{L}(\\{\mathcal{D}_{k}\\}_{k=1}^{M};\theta,\\{H_{k}\\}_{k=1}^{M})=\frac{1}{M}\sum_{k=1}^{M}\mathcal{L}_{k}(\mathcal{D}_{k};\theta,H_{k}),$

where $\mathcal{L}_{k}$ denotes the common loss function in PINNs.
Conventionally, the data for $\mathcal{T}_{k}$ is expressed as
$\mathcal{D}_{k}=\\{\mathcal{D}_{k}^{f},\mathcal{D}_{k}^{b},\mathcal{D}_{k}^{u}\\}$,
where $\mathcal{D}_{k}^{f}=\\{x_{k}^{i},f_{k}^{i}\\}_{i=1}^{N_{k}^{f}}$,
$\mathcal{D}_{k}^{b}=\\{x_{k}^{i},b_{k}^{i}\\}_{i=1}^{N_{k}^{b}}$ and
$\mathcal{D}_{k}^{u}=\\{x_{k}^{i},u_{k}^{i}\\}_{i=1}^{N_{k}^{u}}$, and
$\mathcal{L}_{k}$ as follows:

(3) $\displaystyle\mathcal{L}_{k}(\mathcal{D}_{k};\theta,H_{k})=$
$\displaystyle\frac{w_{k}^{f}}{N_{k}^{f}}\sum_{i=1}^{N_{k}^{f}}||\mathcal{F}_{k}(\hat{u}_{k}(x_{k}^{i}))-f_{k}^{i}||^{2}+\frac{w_{k}^{b}}{N_{k}^{b}}\sum_{i=1}^{N_{k}^{b}}||\mathcal{B}(\hat{u}_{k}(x_{k}^{i}))-b_{k}^{i}||^{2}$
$\displaystyle+$
$\displaystyle\frac{w_{k}^{u}}{N_{k}^{u}}\sum_{i=1}^{N_{k}^{u}}||\hat{u}_{k}(x_{k}^{i})-u_{k}^{i}||^{2}+\mathcal{R}(\theta,H_{k}),$

where $||\cdot||$ represents a properly chosen norm, $\mathcal{R}(\cdot)$ is a
regularization method over the parameters of NNs,
$N_{k}^{f},N_{k}^{b},N_{k}^{u}$ are the numbers of data points for
$f_{k},b_{k},u_{k}$, and $w_{k}^{f},w_{k}^{b},w_{k}^{u}$ are weights to
balance different terms in the loss function.

###  2.2 Generative modeling and normalizing flows (NFs)

As mentioned earlier, MH-PINNs connect $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ by
making two assumptions: (1) the solutions $u_{k},k=1,...,M$ share the same set
of basis functions, $\Phi$; and (2) the corresponding coefficients are samples
of the same random variable, $H$. In [7], $\Phi$ was used as a carrier of
prior knowledge from $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ in downstream physics-
informed learning tasks. In this work, we extend it by utilizing the
information from the head as well by estimating the PDF and a generator of $H$
from its samples, $\\{H_{k}\\}_{k=1}^{M}$, using normalizing flows (NFs). The
interested readers are directed to [32, 22] for reviews of NFs as well as [9,
33, 20] for developments of some popular NFs.

We choose NFs over other commonly used generative models, e.g., generative
adversarial networks (GANs) [13], variational auto-encoders (VAEs) [21], or
diffusion models [16], because the NF serves as both a density estimator and a
generator. The former is able to provide proper regularization in the
downstream few-shot physics-informed learning tasks, while the latter leads to
a physics-informed generative method for approximating stochastic processes.
It is worth noting that in previous works on physics-informed generative
methods [44, 46, 15], NNs are trained by measurements over $u_{k}$, and/or
$f_{k}$, and/or $b_{k}$. Our model, on the other hand, learns through samples
of the head, which is obtained from MTL in the first step. This learning
strategy brings two substantial advantages: (1) flexibility in dealing with
unstructured data, e.g., inconsistent measurements across tasks; (2)
simplicity and controlability of the training by decoupling the physics-
informed learning and the generative modeling.

###  2.3 Prior knowledge utilized in the downstream tasks

Here, we describe details on how to utilize the prior knowledge stored in MH-
PINNs, for downstream few-shot physics-informed learning task,
$\tilde{\mathcal{T}}$, which is defined the same as all other tasks in the
upstream training, but with much fewer measurements. Training of MH-PINNs and
NFs yield a body, $\Phi$, samples of heads, $\\{H_{k}\\}_{k=1}^{M}$, and an
estimated PDF of the head, $\hat{p}(H)\approx p(H)$. In solving
$\tilde{\mathcal{T}}$ with $\tilde{\mathcal{D}}$, we fix the body $\Phi$ and
find the head $\tilde{H}$ that best explains the data $\tilde{\mathcal{D}}$
and the physics in Eq. (1). Noiseless and noisy data are considered in this
paper: for noiseless data, regular NN training is performed on the head for
new tasks to provide deterministic predictions, where the learned PDF of the
head, $\hat{p}(H)$, acts as a regularization term in the loss function; for
noisy data, Bayesian inference is performed on the head as well, in which
$\hat{p}(H)$ denotes the prior distribution. Details are presented in the
following.

####  2.3.1 Regularization in optimization

Limited data in few-shot learning often leads to over-fitting and/or poor
inter-/extrapolation performance. In this regard, regularizing the head
according to its PDF is able to prevent over-fitting and provide additional
prior knowledge for better inter-/extrapolation performance. The optimization
problem is cast as

(4)
$\displaystyle\tilde{H}^{*}=\operatorname*{arg\,min}_{\tilde{H}}\mathcal{L}^{*}(\tilde{\mathcal{D}};\tilde{H}),\text{where
}\mathcal{L}^{*}(\tilde{\mathcal{D}};\tilde{H})=\mathcal{L}(\tilde{\mathcal{D}};\tilde{H})-\alpha\log
p(\tilde{H})$
$\displaystyle\approx\mathcal{L}(\tilde{\mathcal{D}};\tilde{H})-\alpha\log\hat{p}(\tilde{H}),$

where $\mathcal{L}$ is the regular loss function in physics-informed learning
for data $\tilde{\mathcal{D}}$ and parameter $\tilde{H}$, and $\alpha\geq 0$
is the coefficient to adjust the regularization effect. Problem (4) in this
work is solved with gradient descent.

####  2.3.2 Prior distribution in Bayesian inference

As opposed to point estimate obtained by solving the optimization problem (4),
the posterior distribution of the head in $\tilde{\mathcal{T}}$ is obtained
using Bayesian inference. Similar as in [30, 47], the posterior distribution
of $\tilde{H}$ is established as follows:

(5) $p(\tilde{H}|\tilde{\mathcal{D}})\propto
p(\tilde{\mathcal{D}}|\tilde{H})p(\tilde{H})\approx
p(\tilde{\mathcal{D}}|\tilde{H})\hat{p}(\tilde{H}),$

where $p(\tilde{H}|\tilde{\mathcal{D}})$ is the posterior distribution,
$p(\tilde{\mathcal{D}}|\tilde{H})$ is the likelihood distribution, which is
often assumed to be independent Gaussian over all measurements in
$\tilde{\mathcal{D}}$, and $\hat{p}$ is the estimated PDF of the head via NFs.
Intractability of distribution (5) requires approximation methods, among which
Markov chain Monte Carlo methods, such as Hamiltonian Monte Carlo (HMC) [31],
generally provide the most accurate estimation. Moreover, the relatively low
dimension of $\tilde{H}$ also enables the use of Laplace’s approximation (LA)
[18], which is employed in this paper as an alternative to HMC.

##  3 Related works

Deep NNs in recent years have been extensively investigated for solutions of
ODEs/PDEs, SDEs as well as operator learning. Although not explicitly
introduced as MH-PINNs, MH-NNs were first used to solve ODEs/PDEs by [7], in
which MH-PINNs were pre-trained on multiple similar tasks, and then the heads
were discarded while the body was kept and transferred to solving new tasks,
by either least square estimate for linear ODEs/PDEs, or fine-tuning with
gradient descent for nonlinear ones. In [7], a one-shot transfer learning
algorithm for linear problems was proposed but other potential uses of MH-NNs,
e.g., MTL and generative modeling, were not discussed, as opposed to the work
presented herein. Furthermore, [7] focused only on fast and deterministic
predictions with high accuracy using sufficient clean data, while in this
paper, we study the applicability of MH-NNs to few-shot physics-informed
learning as well, where data is insufficient and/or noisy, and address such
cases with UQ. We note that MH-NN was also used as a multi-output NN in [45],
which, however, focused on solving single tasks and obtaining uncertainties.

Generative modeling in the context of scientific computing has also been
studied recently, and a few attempts for adopting deep generative NNs to SciML
problems have been made in [44, 46, 15], most of which have focused on
approximating stochastic processes and on solving SDEs. We propose a new
physics-informed generative method, as an alternative to the current ones,
using MH-PINNs, and test it in approximating stochastic processes. In this
regard, our method is functionally the same as the current ones, but
technically different. All previous methods address physics-informed
generative modeling using end-to-end learning strategies by coupling two
dissimilar types of learning, physics-informed learning and generative
modeling, which may be problematic for implementation and usage in practice
when either type of learning becomes more complicated. Our method, on the
other hand, addresses the problem in an entirely new angle by decoupling those
two: physics-informed learning is performed first and is followed by learning
generators. To this end, our method is a two-step method, and with the help of
well-developed algorithms from both fields, our method has advantages both in
flexibility and simplicity in implementation.

##  4 Results

In this section, we test our method using five benchmarks. The first one is a
pedagogical function regression, in which we aim to demonstrate the basic
applicability and capabilities of our method, showing the importance of
incorporating the distribution of the head in the downstream tasks and in
obtaining results with or without uncertainty. The second example is a
nonlinear ODE system, in which we test our method in approximating stochastic
processes through a differential operator, compare different NFs, and
eventually compare our method with another well-known physics-informed
generative model, physics-informed GANs (PI-GANs) [44] in generative modeling.
The third is a 1-D nonlinear reaction-diffusion equation, the fourth is a 2-D
nonlinear Allen-Cahn equation, and the fifth is the 2-D stochastic Helmholtz
equation with $20$ dimensions. In all examples unless stated otherwise, data
for $\\{\mathcal{D}_{k}\\}_{k=1}^{M}$ are noise-free and task-wisely
sufficient, while $\tilde{\mathcal{D}}$ in downstream tasks is insufficient,
which makes the downstream tasks of the few-shot type. In addition, except for
the first example, results from Bayesian inference are obtained by employing
HMC, and the predicted mean denoted as $\mu$ and predicted standard deviation
denoted as $\sigma$ are computed from the posterior samples of functions or
unknown parameters. The predicted uncertainty is defined as $2\sigma$ in this
paper.

###  4.1 Function approximation

We start with a function regression problem using only data and no physics,
which is a degenerate instance of Eq. (1) with $\mathcal{F}_{k}$ being fixed
as an identity operator, no $\mathcal{B}$ and $b_{k}$, and $u_{k}=f_{k}$ being
task-specific. In this case, $\mathcal{D}_{k}$ and $\tilde{\mathcal{D}}$ are
given as $\\{(x_{k}^{i},f_{k}^{i})\\}_{i=1}^{N_{k}}$ and
$\\{(x^{i},f^{i})\\}_{i=1}^{N}$, respectively, and $\mathcal{T}_{k}$ and
$\tilde{\mathcal{T}}$ are defined as approximating functions $f_{k}$ and
$\tilde{f}$ from $\mathcal{D}_{k}$ and $\tilde{\mathcal{D}}$, respectively.
The stochastic function $f$ in this example is defined as follows:

(6) $\displaystyle f(x)=A\cos(\omega x)+2\beta x,x\in[-1,1],$ $\displaystyle
A\sim U[1,3),\omega\sim U[2\pi,4\pi),\mathcal{P}(\beta=\pm 1)=0.5,$

where $U$ stands for uniform distribution and $\mathcal{P}(\Xi)$ is defined as
the probability of the event $\Xi$. Our goal is to approximate $f$ from data
with MH-NNs and NFs, and solving the downstream few-shot regression tasks
$\tilde{\mathcal{T}}$ as well, in which two functions, $2\cos(2\pi x)+2x$ and
$2\cos(4\pi x)-2x$, are regressed from $4$ and $5$ measurements equidistantly
distributed on $[-0.9,-0.1]$, respectively.

For the training of MH-NNs and NFs, $1,000$ $f$ subject to Eq. (6) are
sampled, each of which forms a regression task with $40$ measurements sampled
equidistantly on $[-1,1]$ as data. Samples of $f$ for training are displayed
in Fig. 3(a). Both noiseless and noisy data cases are considered in the few-
shot regression tasks. As described in Sec. 2.3, the former is solved by fine-
tuning the head using gradient descent, while the latter is solved by
estimating the posterior distribution (Eq. (5)) using HMC and LA. The noise
$\varepsilon$ is assumed to be independent additive Gaussian noise with scale
$0.2$, i.e., $\varepsilon\sim N(0,0.2^{2})$. In the downstream few-shot
regression tasks we compare our method with two other approaches, the transfer
learning (TL) method from [7], which only transfers the body, $\Phi$, and the
regular NN method, in which no prior knowledge is employed and all parameters
of NN are trained.

Figure 3: Results for approximating the stochastic function defined in Eq. (6)
and solving the downstream few-shot regression tasks. (a) Left: $1,000$
samples generated from the exact distribution; middle: $1,000$ samples
generated from the learned generator; right: statistics computed from samples,
in which we refer to the interval of mean $\pm$ 2 standard deviations as
bound. (b)/(c) Results for the downstream tasks. Left: results for noiseless
cases using our method, the transfer learning (TL) method in [7], and regular
NN method; middle: results for noisy case using our method with HMC for
posterior estimate; right: results for the same noisy case using our method
with LA for posterior estimate.

Results for approximating $f$ and solving the downstream tasks are presented
in Fig. 3. As shown in Fig. 3(a), our method approximates the stochastic
function $f$ well, demonstrating the capability of MH-NNs in generative
modeling. In solving the downstream tasks with noiseless data, $L_{2}$
regularization is imposed in the TL method and regular NN method, to prevent
over-fitting when only $4$ or $5$ measurements are available. As we can see
from Figs. 3(b) and (c), our approach yields accurate predictions and performs
significantly better than the other two in both tasks, particularly in the
region where there are no measurements. By comparing our approach with the NN
method, we can see that prior knowledge of $f$ is learned from
$\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ and transferred successfully to new
downstream tasks. By comparing our approach with the TL method, we can see
that the prior knowledge is stored in both the body and (the distribution of)
the head. For the noisy cases, it is shown that, for both tasks and both
posterior estimating methods, the predictions are accurate and trustworthy:
the predicted means agree with the references and the errors are bounded by
the predicted uncertainties. It is worth noting that the predicted
uncertainties do not develop in the interval $[0,1]$ and show periodic
patterns, even if there are no measurements. That is because an informative
prior, which is learned by MH-NNs and NFs, is imposed on the head in Bayesian
inference.

The target functions in downstream tasks considered previously are chosen to
be in-distribution. They are regressed well with insufficient data, mainly
because they belong to the space of functions, on which the generator is
trained. However, when functions in the downstream tasks are out-of-
distribution (OOD), our approach fails to produce good predictions, even if
the data is sufficient, as shown in Fig. 4. Here, the target function is
chosen to be $2\cos(4.5\pi)+x$ with both $\omega$ and $\beta$ being OOD.
Fluctuations are predicted but do not match the reference. In Fig. 4, we can
further see that when data is sufficient, a NN trained from scratch
significantly outperforms our approach, showing that, for OOD functions, the
more we rely on the learned regularization, which is indicated by the value of
$\alpha$ in Eq. (4), the more erroneous the prediction is.

Figure 4: Results for regression on an out-of-distribution function. Left:
few-shot regression with clean data using our approach; middle: few-shot
regression with noisy data using our approach with HMC for posterior estimate;
right: regression with sufficient clean data using regular NN method and our
approach with different regularization terms, $\alpha$ in Eq. (4).

###  4.2 Nonlinear ODE system

In this example, we consider the following ODE system [28], which describes
the motion of a pendulum with an external force:

(7) $\displaystyle\frac{du_{1}}{dt}$ $\displaystyle=u_{2},$
$\displaystyle\frac{du_{2}}{dt}$ $\displaystyle=-\lambda\sin(u_{1})+f(t),$

with initial condition $u_{1}(0)=u_{2}(0)=0$. In Eq. (7), $f$ is the external
force and $\lambda$ is a constant. Here, to demonstrate and study the
capability of our method in generative modeling, we first consider the case
where $f$ is a Gaussian process and $\lambda=1$ is known, which is referred to
as the forward problem. Different from previous studies [44, 46, 15], in which
a stochastic process is approximated directly by the output of NNs, in this
example we place the differential operator right after NNs and approximate the
stochastic process $f$ as the source term in Eq. (1). We also test our method
on the inverse problem, where the values of $\lambda$ in Eq. (7) are unknown
in $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ and $\tilde{\mathcal{T}}$. The forward
problem corresponds to Eq. (1) with $\mathcal{F}_{k},b_{k}$ being the same for
all tasks and $u_{k},f_{k}$ being task-specific, while the inverse problem
corresponds to Eq. (1) with $b_{k}$ being the same and $u_{k},f_{k}$, and the
differential operator $F_{k}$ being different as a consequence of task-
specific $\lambda$.

####  4.2.1 Forward problem

We first assume $\lambda=1$ in Eq. (7) is known, and the data on $f$ is
available, i.e., $\mathcal{D}_{k}=\\{(x_{k}^{i},f_{k}^{i})\\}_{i=1}^{N_{k}}$,
$k=1,...,M$. As described before, we employ MH-PINNs to solve
$\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ all at once and then employ NFs to learn the
distribution of the head. Consequently, we obtain generators of $f$ and $u$.
In this case, $f$ is assumed to be a Gaussian process with squared kernel
function:

(8) $\displaystyle f(t)\sim
GP(0,K),t\in[0,1],K(x,x^{\prime})=\exp(-\frac{|x-x^{\prime}|^{2}}{2l^{2}}),$

where the correlation length $l$ is set to $0.1,0.25,0.5$. As discussed in
Sec. 2.2, many types of NFs have been developed in the past decade for
generative modeling and density estimate. MH-PINNs, when used as generators,
are compatible with all NFs. In this regard, we compare three popular NFs,
RealNVP [9], MAF [33] and IAF [20], and eventually compare MH-PINNs (with NFs)
against PI-GAN [44] in approximating the Gaussian process defined in Eq. (8)
with different correlation lengths.

For the training of MH-PINNs and NFs as well as PI-GANs, $2,000$ $f$ are
sampled with respect to Eq. (8), each of which forms a physics-informed
regression task with $65$ measurements of $f$ equidistantly sampled on $[0,1]$
as data. Notice that the ODE system in Eq. (7) can be rewritten in a simpler
format as follows:

(9) $\displaystyle u_{tt}=-\lambda\sin(u)+f(t),t\in[0,1],$

with initial conditions $u(0)=u_{t}(0)=0$. Hence, we choose to use Eq. (9) to
build the loss function for physics-informed learning.

Results for comparisons are shown in Fig. 5 and Table 1. From the spectral
analysis of the approximated Gaussian processes shown in Fig. 5, we can see
that MH-PINNs with MAF and IAF are comparable with PI-GANs while MH-PINNs with
RealNVP fall marginally behind. As shown in Table 1, the computational costs
of MH-PINNs with MAF and RealNVP are significantly lower than PI-GANs, while
MH-PINNs with IAF is more expensive than PI-GANs. We note that in this example
we also record the computational cost for sampling using different generators.
As shown in Table 1, PI-GANs are significantly faster in generating samples.
That is because, generally, GANs require relatively shallow NNs as opposite to
NFs, for which a deep architecture is needed to keep up the expressivity.
Among three NFs, IAF is the fastest in sampling while MAF is the lowest, as
opposed to training, which is consistent with the properties of those two NFs:
MAF is slow for the forward pass, which is used to generate samples, and fast
for the inverse pass, which is used to compute the density, while IAF is the
opposite. Despite the fact that MAF is slow in sampling, considering its fast
training and good performance, we equip MH-PINNs with MAF as the density
estimator and the generator for all other examples in this paper.

Figure 5: Approximating Gaussian processes as the source term in Eq. (7) using different models: spectra of the correlation structure for the learned generators, for different correlation lengths, $l$. The covariance matrix is constructed using $10,000$ generated samples, and eigen-values are averaged over 10 generators trained independently. | MAF | IAF | RealNVP | PI-GAN
---|---|---|---|---
Phase 1 |  $134$s |  $134$s |  $134$s | N/A
Phase 2 |  $252$s |  $3939$s |  $245$s | N/A
Total |  $386$s |  $4073$s |  $379$s |  $3243$s
Sampling |  $1.98\times 10^{-1}$s |  $1.48\times 10^{-2}$s |  $1.50\times 10^{-2}$s |  $2.29\times 10^{-3}$s
Table 1: Computational time for different models to approximate Gaussian
process with correlation length $l=0.1$. The MH-PINN method is a two-step
method and hence its computation is decomposed into two parts: training MH-
PINNs referred to as phase 1 and training NFs referred to as phase 2. Sampling
time is defined to be the average time needed to generate $10,000$ samples of
$u$.

####  4.2.2 Inverse problem

Next, we assume $\lambda$ in Eq. (7) is unknown, and some measurements of $u$
are available, in addition to $f$, i.e.,
$\mathcal{D}_{k}=\\{\\{x_{k}^{i},f_{k}^{i}\\}_{i=1}^{N_{k}^{f}},\\{x_{k}^{i},u_{k}^{i}\\}_{i=1}^{N_{k}^{u}}\\}$.
MH-PINNs are first employed to infer $u_{k}$ as well as $\lambda_{k}$ from
data $\mathcal{D}_{k}$ and physics, and NFs are employed afterwards to learn
from samples of $H_{k}$ and $\lambda_{k}$. To this end, the generative model
is for the joint distribution of $u,f$ and $\lambda$. Here, we assume $f$
follows a truncated Karhuen-Loeve (KL)-expansion, with $5$ leading terms, of
the Gaussian process with squared kernel function and correlation length
$0.1$, and for each task $\mathcal{T}_{k}$,
$\lambda_{k}=\frac{1}{2}\exp(\int_{[0,1]}f^{2}_{k}(t)dt)$. As for the
downstream task $\tilde{\mathcal{T}}$, the target is to infer $u$ and
$\lambda$ from insufficient data of $u$ and $f$.

For the training of MH-PINNs and NFs, $2,000$ samples of $f$ are generated and
displayed in Fig. 6(a). For each task, we assume $33$ measurements of $f_{k}$
and $9$ measurements of $u_{k}$, equidistantly distributed on $[0,1]$, are
available, and initial conditions are hard-encoded in NN modeling. For the
downstream task, we assume $1$ random measurement of $u$ and $8$ random
measurements of $f$ are available with hard-encoded initial conditions as
well. For the case with noisy measurements, we assume the noises
$\varepsilon_{f}$ and $\varepsilon_{u}$ to be additive Gaussian, with $0.05$
noise scale for measurements of $f$ and $0.005$ noise scale for measurements
of $u$, respectively, i.e. $\varepsilon_{f}\sim N(0,0.05^{2})$ and
$\varepsilon_{u}\sim N(0,0.005^{2})$. The reference solution as well as the
clean data of $u_{k}$ are generated by solving Eq. (7) for each task
$\mathcal{T}_{k}$, with corresponding $f_{k}$ and $\lambda_{k}$ using Matlab
ode45.

Results are shown in Fig. 6 and Table 2, from which we can see our method is
able to approximate the stochastic process as a source term well and produce
accurate and trustworthy predictions, for $u,f$ and also $\lambda$ in the
downstream task with limited data, in both noiseless and noisy cases. As
shown, the PINN method yields unacceptable estimate over both $u$ and
$\lambda$ due to lack of data, while our approach is of much higher accuracy
by integrating prior knowledge from $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ with MH-
PINNs.

Figure 6: Results for the inverse problem of the ODE system (7), with initial conditions hard-encoded in NN modeling. (a) Left: $1,000$ samples of $f$, generated from the exact distribution; middle: $1,000$ samples of $f$, generated from the learned generator; right: statistics computed from samples. The bound is defined as the same as in the caption of Fig. 3. (b) Predicted $f$ and $u$ using PINNs and our approach, for the downstream inverse problem with noiseless data. (c) Predicted $f$, $u$ and $\lambda$ with uncertainties using our approach, for the downstream inverse problem with noisy data. The predicted mean and standard deviation of $\lambda$ is $2.4663$ and $0.1501$, while the reference value is $2.3609$. | PINN | MH-PINN
---|---|---
$\lambda$ | $0.8440$ | $2.5428$
Error (%) | $63.99$ | $1.21$
Table 2: Estimate of $\lambda$ and $L_{2}$ relative error of $u$ for the
downstream inverse problem on Eq. (7) with clean data, using our approach and
the regular PINN method. The reference value for $\lambda$ is $2.3609$.

###  4.3 1-D nonlinear reaction-diffusion equation

We now test our method on a 1-D nonlinear time-dependent reaction-diffusion
equation, which is commonly referred to as Fisher’s equation [2]:

(10) $\displaystyle u_{t}=Du_{xx}+ku(1-u),t\in[0,1],x\in[-1,1],$ (11)
$\displaystyle u(t,-1)=u(t,1)=0,t\in[0,1],$ (12) $\displaystyle
u(0,x)=u_{0}(x),x\in[-1,1],$

where $D=0.1,k=0.1$ and $u_{0}(x)$ is the initial condition function. In this
example, we assume that the initial condition function is a stochastic process
with the following distribution:

(13) $\displaystyle
u_{0}(x)=\frac{(x^{2}-1)}{5}\sum_{j=1}^{5}\xi_{j}(\cos^{2}(jx)-1),x\in[-1,1],$

where $\xi_{j},j=1,...,5$ are independent and identically distributed (i.i.d.)
random variables subject to uniform distribution on $[0,1)$, i.e.,
$\xi_{j}\sim U[0,1)$. Unlike previous examples, the stochasticity comes from
the initial condition rather than the source term. This example corresponds to
Eq. (1) with $\mathcal{F}_{k}$ and $f_{k}$ being the same for all tasks, and
$u_{k}$ and $b_{k}$ being task-specific. In addition to measurements of
$u_{0}$, points on which the PDE residuals are computed are also required in
both $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ and $\tilde{\mathcal{T}}$. Hence, the
data is
$\mathcal{D}_{k}=\\{\\{(t_{k}^{i},x_{k}^{i}),0\\}_{i=1}^{N_{k}^{f}},\\{(0,x_{k}^{i}),b_{k}^{i}\\}_{i=1}^{N_{k}^{b}}\\}.$

For the training, $2,000$ samples of $u_{0}(x)$ are generated, displayed in
Fig. 7(a), and each sample forms a physics-informed regression task with $41$
measurements of $u_{0}$ equidistantly sampled on $[-1,1]$ as data for initial
condition. Besides, for all tasks, a uniform mesh $21\times 41$ on temporal-
spatial domain $[0,1]\times[-1,1]$ is used to compute the PDE residual loss.
For the downstream tasks $\tilde{\mathcal{T}}$, $5$ random measurements of
$u_{0}$ are available and the same uniform mesh is applied. The boundary
conditions are hard-encoded in NN modeling in both
$\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ and $\tilde{\mathcal{T}}$. For the noisy
case, the noise $\varepsilon$ is assumed to be independent additive Gaussian
noise with $0.02$ noise scale for both measurements of $u_{0}$ and the PDE
residual, i.e., $\varepsilon\sim N(0,0.02^{2})$. Results are presented in Fig.
7 and Table 3. We can see that our method estimates a good generator of the
stochastic processes from data and physics, which provides informative prior
knowledge in the downstream few-shot physics-informed regression tasks. The
prediction is accurate in both noiseless and noisy cases, and the errors in
the noisy case are bounded by the predicted uncertainty. The $L_{2}$ error of
$u$, shown in Table 3, indicates that our approach outperforms the PINN method
by a significant amount, hence demonstrating the effectiveness of bringing
prior knowledge into solving similar tasks.

Figure 7: Generator learning and few-shot physics-informed learning on 1-D time-dependent reaction-diffusion equation (10), with boundary conditions hard-encoded in NN modeling. (a) Left: $1,000$ training samples of $u_{0}$; middle: $1,000$ samples of $u(0,\cdot)$ from the learned generator; right: statistics computed from samples. The bound is defined as the same as in the caption of Fig. 3. (b) Predicted $u$ at $t=0,0.5,1$ using our approach and the PINN method with noiseless measurements. (c) Predicted mean and uncertainty of $u$ at $t=0,0.5,1$ using our approach with HMC for posterior estimate, with noisy measurements. | PINN | MH-PINN
---|---|---
Error (%) | $78.77$ | $0.22$
Table 3: $L_{2}$ relative error of $u$ for the downstream few-shot physics-
informed learning task on Eq. (10) with clean data of $u_{0}$ using our
approach and the PINN method.

###  4.4 2-D nonlinear Allen-Cahn equation

We now move to a 2-D steady nonlinear Allen-Cahn equation with Dirichlet
boundary conditions [43]:

(14) $\displaystyle\lambda\Delta u+u(u^{2}-1)$ $\displaystyle=f,x,y\in[0,1],$
(15) $\displaystyle u(x,0)=u(x,1)$ $\displaystyle=u(0,y)=u(1,y)=0,$

where $\lambda=0.1$ is a constant and $f$ is the source term. Here, we impose
a distribution to $f$, which is derived from Eq. (14) and the following
distribution of the solution $u$:

(16) $\displaystyle u(x,y)=\frac{1}{5}\sum_{j=1}^{5}\xi_{j}\frac{\sin(j\pi
x)\sin(j\pi y)}{j^{2}\pi^{2}},x,y\in[0,1],$

with i.i.d. random variables $\xi_{j},j=1,...,5,$ subject to uniform
distribution, i.e. $\xi_{j}\sim U[0,1)$. In this example, we wish to use our
method to learn generators of both $u$ and $f$ from data of $f$ and physics in
Eq. (14), and use it to solve the downstream task $\tilde{\mathcal{T}}$ with
insufficient data $\tilde{\mathcal{D}}$. This example corresponds to Eq. (1)
with $\mathcal{F}_{k},b_{k}$ being the same among tasks and $f_{k},u_{k}$
being task-specific, and the data is
$\mathcal{D}_{k}=\\{(x_{k}^{i},y_{k}^{i}),f_{k}^{i}\\}_{i=1}^{N_{k}^{f}}$.

To train the MH-PINNs and NFs, we sample $2,000$ $f$ from its distribution,
each of which is resolved with a $51\times 51$ uniform mesh on 2-D spatial
domain $[0,1]\times[0,1]$. As for the downstream task, $100$ random
measurements of $f$ on the uniform mesh are assumed to be available. The noise
is assumed to be independent additive Gaussian noise with $0.05$ noise scale.
In both $\mathcal{T}_{k}$ and $\tilde{\mathcal{T}}$, the boundary conditions
are hard-encoded in NN modeling. Results as well as the locations of the
measurements are presented in Fig. 8 and Table 4. Similar to all previous
examples, our approach delivers accurate and trustworthy predictions, showing
that prior knowledge is learned and transferred well in both deterministic and
Bayesian inferences.

Figure 8: Results for few-shot physics-informed learning on the 2-D nonlinear Allen-Cahn equation Eq. (14) with noisy measurements of $f$. Predicted mean $\mu$ and standard deviation $\sigma$ are computed over $1,000$ posterior samples from HMC. The absolute error is defined as the absolute value of difference between the reference and $\mu$. Black crosses represent the locations of the measurements on $f$. | PINN | MH-PINN
---|---|---
Error (%) | $12.82$ | $0.30$
Table 4: $L_{2}$ relative error of $u$ for the downstream few-shot physics-
informed learning task on Eq. (14) with clean data of $f$, using our approach
and the PINN method.

###  4.5 2-D Stochastic Helmholtz equation

The last example we test in this paper is the 2-D Helmholtz equation with
stochastic source term and Dirichlet boundary conditions [35]:

(17) $\displaystyle(\lambda^{2}-\nabla^{2})u$
$\displaystyle=f,x,y\in[0,2\pi],$ (18) $\displaystyle u(x,0)=u(x,2\pi)$
$\displaystyle=u(0,y)=u(2\pi,y)=0,$

where $\lambda^{2}$ is the Helmholtz constant and $f$ is defined as follows:

(19) $\displaystyle
f(x,y)=\frac{2}{d}\\{\sum_{i=1}^{d/4}\xi_{i}\sin(ix)+\xi_{i+d}\cos(ix)+\xi_{i+2d}\sin(iy)+\xi_{i+3d}\cos(iy)\\},$

where $\xi_{j},j=1,...,d$ are i.i.d. random variables subject to uniform
distribution $U[0,1)$ and $d$ represents the dimension of the randomness. For
demonstration purposes, we consider the case where $d=20$ in this paper,
unlike the one in [35] with $d=100$.

The first case we study is the forward problem with $\lambda^{2}=1$ known.
This setup corresponds to Eq. (1) with $\mathcal{F}_{k},b_{k}$ being shared
among tasks and $u_{k},f_{k}$ being task-specific. Next, we study the inverse
problem with unknown $\lambda$, where data on $u$ and $f$ are available, which
corresponds to Eq. (1) with only $b_{k}$ being the same and $u_{k},f_{k}$ and
operator $\mathcal{F}_{k}$ being task-specific. The downstream tasks are
defined as the same as $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$ in both cases, but
with fewer measurements.

For both the forward and inverse problems, $10,000$ $f$ are sampled from its
distribution, and hence $10,000$ tasks are solved with MH-PINNs with boundary
conditions hard-encoded in NN modeling. We display the samples of a slice of
$f$ in Fig. 9(a). For the forward problem, $\mathcal{D}_{k}$ only contains
measurements of the source term $f_{k}$, i.e.
$\mathcal{D}_{k}=\\{\\{(x_{k}^{i},y_{k}^{i}),f_{k}^{i}\\}_{i=1}^{N_{k}^{f}}\\}$,
while for the inverse problem $\mathcal{D}_{k}$ also contains measurements of
the sought solution $u_{k}$:
$\mathcal{D}_{k}=\\{\\{(x_{k}^{i},y_{k}^{i}),f_{k}^{i}\\}_{i=1}^{N_{k}^{f}},\\{(x_{k}^{i},y_{k}^{i}),u_{k}^{i}\\}_{i=1}^{N_{k}^{u}}\\}$.
For the training in the forward problem, each sample of $f$ is resolved by a
$50\times 50$ uniform mesh on 2-D spatial domain $(0,2\pi)\times(0,2\pi)$ with
boundary excluded. For the inverse problem, the same $10,000$ samples of $f$
are used, but this time they are resolved with a $21\times 21$ uniform mesh.
In addition, for each task $\mathcal{T}_{k}$, measurements of $u_{k}$ on a
$6\times 6$ uniform mesh are available. The reference solution and
measurements of $u$ are generated by solving Eq. (17) with
$\lambda_{k}^{2}=\int_{[0,2\pi]^{2}}f_{k}^{2}(x,y)dxdy$ using the finite
difference method with five-point stencil. For the downstream tasks, $100$
random measurements of $f$ are available for the forward problem, and $50$
random measurements of $f$ and $10$ random measurements of $u$ are available
for the inverse problem. The noise is assumed to be independent additive
Gaussian noise with $0.05$ noise scale.

Results are displayed in Tables 5 and 6, and Figs. 9 and 10. As shown, the
learned generator is able to produce samples of $f$ with high quality as well
as providing informative prior knowledge for the downstream tasks, in both the
forward and inverse problems. As for the noisy case with Bayesian inference
and UQ, the predicted means agree with the references and the absolute errors
are mostly bounded by the predicted uncertainties. The effectiveness of our
approach for few-shot physics-informed learning and the applicability to both
deterministic optimization and Bayesian inference have been consistently
demonstrated in the past five examples.

Figure 9: Generator learning and few-shot physics-informed learning on the stochastic Helmholtz equation (17). (a) Left: $1,000$ training samples of a slice of $f$ at $y=\pi$; middle: $1,000$ samples of a slice of $f$ at $y=\pi$ from the learned generator; right: statistics computed from samples. (b)/(c) Results for the downstream forward problem with $100$ random noisy measurements on $f$, using our approach with HMC. From left to right are reference, predicted mean, absolute error, and predicted uncertainty of $f$/$u$. Black crosses represent the locations of the measurements of $f$. | PINN | MH-PINN
---|---|---
Error (%) | $21.14$ | $1.12$
Table 5: $L_{2}$ relative error of $u$ for the downstream forward problem on
Eq. (17) with clean data of $f$, using our approach and the PINN method.

Figure 10: Results for the downstream inverse problem on the stochastic Helmholtz equation (17), with $50$ random noisy measurements of $f$ and $10$ random noisy measurements of $u$. $\lambda$ is estimated as $1.0785\pm 0.0307$ in the format of predicted mean $\pm$ predicted standard deviation, while the reference value is $1.0042$. (a)/(b) From left to right are reference, predicted mean, absolute error, and predicted uncertainty of $f$/$u$. Black crosses represent locations of the measurements of $f$ or $u$. | PINN | MH-PINN
---|---|---
$\lambda$ | $1.9328$ | $1.0170$
Error (%) | $59.92$ | $2.58$
Table 6: Estimate of $\lambda$ and $L_{2}$ relative error of $u$ for the
downstream inverse problem on Eq. (17) with clean data. The reference value of
$\lambda$ is $1.0042$.

##  5 Multi-task learning with multi-head neural networks

So far we have mostly focused on using MH-NNs together with NFs to estimate
stochastic generators and learn informative prior knowledge from
$\\{\mathcal{T}_{k}\\}_{k=1}^{M}$. This was achieved by first training MH-NNs
in a MTL fashion and then training NFs to estimate the PDF of the head.
Intuitively, the capability of MH-NNs when trained in MTL in capturing shared
information is the key to the success in generative modeling and few-shot
learning. For physics-informed MTL with MH-PINNs, ODEs/PDEs are solved
simultaneously, and assuming the solutions to share the same set of basis
functions gives us samples of the set of coefficients, which enables the
generative modeling, followed by few-shot learning, which is the whole point
of the method proposed in this paper. However, the cost and/or the benefit of
imposing the same set of basis functions to all solutions have not been
explicitly discussed yet. On one hand, the shared body relates the training of
tasks, which may be helpful if tasks are similar in certain ways. On the other
hand, forcing all solutions to share the same basis functions may also be
harmful when they behave differently. In particular, for tasks with sufficient
data and physics, forcing them to share the same body with all other tasks may
act as a negative regularization, and single-task learning (STL) may
outperform MTL in terms of the prediction accuracy in those specific tasks. In
this section, we investigate the effect of MTL using MH-NNs and provide
preliminary results and analysis by revisiting the simple function
approximation example in Sec. 4.1, which, hopefully, could provide useful
information and insight for future more rigorous research.

###  5.1 Basis function learning and synergistic learning

As discussed before, the quality and behavior of basis functions learned in
MTL are crucial to generative modeling and learning the relation and the
representative information of tasks $\\{\mathcal{T}_{k}\\}_{k=1}^{M}$. We
consistently noticed from numerical examples that the initialization of the
head in MH-NNs has great impact on the average accuracy of MTL, the learning
of the basis functions, and the distribution of the head. Here, we test three
initialization strategies, random normal method with $0.05$ standard deviation
referred to as RN (0.05), Glorot uniform method [12] referred to as GU, and
random normal method with $1$ standard deviation referred to as RN (1). In the
downstream few-shot learning tasks, we fine-tune the head without the learned
PDF, which is in fact the TL method from [7], by which the information from
the distribution of the head is excluded and the prediction accuracy is fully
determined by the level of prior knowledge contained in the basis functions.

As shown in Fig. 11, RN ($0.05$) yields the least informative basis functions,
whose behavior is dominated by the hyperbolic tangent activation function of
NNs. This is further demonstrated in the downstream few-shot learning tasks
using the TL method. It also provides the worst prediction accuracy on average
in MTL, as presented in Table 8. GN and RN ($1$) perform similarly. Plots of
some basis functions seemingly indicate that RN ($1$) yields better basis
functions, whose behaviors are more similar to the family of functions
displayed in Fig. 11(b), which, however, does not necessarily imply richer
prior knowledge in the downstream tasks, as shown in Fig. 11(c).

It is shown empirically that compared to other two initialization strategies,
MH-NNs with RN ($0.05$) does not deliver accurate MTL nor synergistic learning
in basis functions. However, we noticed that, in generative modeling, it
performs significantly better in terms of accuracy and convergence speed. As
shown in Fig. 11(d), samples from the learned generator are of higher quality.
We consistently found that initializing heads with relatively small values
often led to easy and fast training of NFs and accurate learning of the
generative models. We conjecture that this happens because MH-NNs in MTL tend
to contain the representative and informative information in the heads when
heads are initialized with small values, while contain it in the basis
functions when heads are initialized with relatively large values.

|  RN ($0.05$) | GU |  RN ($1$)
---|---|---|---
Error (%) | $0.8373\pm 0.2341$ | $0.1907\pm 0.0690$ | $0.3131\pm 0.0937$
Table 7: $L_{2}$ relative errors, from MTL, for $1,000$ tasks, using different
initialization methods. The errors are displayed in the format of mean $\pm$
standard deviation, computed over all tasks.

Figure 11: The effect of different initialization methods of the head, in
basis functions learning, few-shot learning, and generator learning. (a)
Samples of $20$ basis functions from MH-NNs, trained for approximating $1,000$
$f$ generated from Eq. (6), using, from left to right, RN ($0.05$), GU and RN
($1$) initialization methods. (b) $1,000$ training samples of $f$. (c) Results
for two downstream few-shot regression tasks, using TL method without
regularization informed by the learned PDF, as opposite to the proposed
approach. (d) Results for generator learning, using, from left to right, RN
($0.05$), GU and RN ($1$) initialization methods.

###  5.2 Multi-Task Learning (MTL) versus Single-Task Learning (STL)

As discussed earlier, MTL with MH-NNs does not necessarily result in
synergistic learning nor higher accuracy for all tasks on average. Here, we
use again the function approximation example in Sec. 4.1, to investigate the
effectiveness of MTL with MH-NNs, as compared to STL. The first case we
consider here is assuming that the data is sufficient. For that, we randomly
choose $100$ out of the $1,000$ training samples, each one of which is
approximated by a NN trained independently, and compare the results with MH-
NNs in terms of prediction accuracy. Note that in this case, a MH-NN is
trained on $1,000$ functions as before and tested on the chosen $100$
functions, while a single-head NN with the same architecture is trained on
$100$ functions directly. Results are shown in Table 8, from which it is
verified empirically that MTL is outperformed by STL under certain
circumstances, e.g., when the random normal initialization methods are used.

|  RN ($0.05$) | GU |  RN ($1$) | STL
---|---|---|---|---
Error (%) | $0.7575\pm 0.2477$ | $0.1362\pm 0.0259$ | $0.3664\pm 0.1031$ | $0.2102\pm 0.0794$
Table 8: $L_{2}$ relative errors of $f$, from MTL with MH-NNs and STL with
NNs, on $100$ tasks. Different initialization methods are used for the heads
in MH-NNs. The errors are displayed in the format of mean $\pm$ standard
deviation, computed over all $100$ tasks.

The second case we consider is assuming that the data is sufficient for some
tasks while insufficient for other tasks. For that, we split equally the
$1,000$ tasks into two subsets of tasks. For the first $500$ tasks, we assume
we only have $10$ measurements randomly sampled on $[-1,1]$, while for the
other $500$ tasks, we assume we have full $40$ measurements equidistantly
distributed on $[-1,1]$. MTL with MH-NNs is performed on those $1,000$
regression tasks all at once, and the tasks are treated as equal. The results
are presented in Table 9 and Fig. 12. We can see that, compared to STL, MTL
improves the prediction accuracy on tasks with insufficient data, providing
empirical evidence of synergistic learning. Also, interestingly, RN ($1$)
initialization method, which yields the worst generative models, performs the
best among all three, which agrees with our previous conjecture on the basis
functions learning with MH-NNs, that heads initialized with large values tend
to force representative and informative information to be encoded in the basis
functions.

|  RN ($0.05$) | GU |  RN ($1$)
---|---|---|---
Error (%) | $63.60\pm 24.08$ | $40.49\pm 20.49$ | $16.91\pm 11.08$
Table 9: $L_{2}$ relative errors of $f$, from MTL with MH-NNs, on $500$ tasks
equipped with insufficient data. The errors are displayed in the format of
mean $\pm$ standard deviation, computed over all $500$ tasks.

Figure 12: Results for 3 tasks with insufficient data from MTL with MH-NNs,
using different initialization methods over the head, and from STL with NNs
with the same architecture. We note that tasks with sufficient data and tasks
with insufficient data are treated equally in MTL.

##  6 Discussion

We have developed multi-head neural networks (MH-NNs) for physics-informed
machine learning, and proposed multi-head physics-informed neural networks
(MH-PINNs) as a new method, implemented in the L-HYDRA code. The primary focus
of this work is on MH-NNs and MH-PINNs for various learning problems in
scientific machine learning, including multi-task learning (MTL), stochastic
processes approximation, and few-shot regression learning. We first formulated
the problem in Eq. (1), introduced the architecture design of MH-PINNs, and
proposed a method to transform MH-NNs and MH-PINNs to generative models with
the help of normalizing flows (NFs) for density estimation and generative
modeling. We then studied the applicability and capabilities of MH-PINNs in
solving ordinary/paritial differential equations (ODEs/PDEs) as well as
approximating stochastic processes. We completed the paper with preliminary
and empirical explorations of MH-NNs in synergistic learning, and examined the
potential benefits and cost of MTL with MH-NNs.

This paper can be used in various ways: it proposes a NN approach for MTL in
solving ODEs/PDEs; it provides a new approach to approximate stochastic
processes; it presents a method to address few-shot physics-informed learning
problems, which are often encountered in the context of meta-learning and
transfer learning; it contains a systematic study of applying MH-NNs to
scientific computing problems; it presents the first empirical evidence of
synergistic learning.

However, there are a few major problems on MH-NNs we did not address, one of
which is the expressivity of MH-NNs, or more generally hard-parameter sharing
NNs in approximating complicated stochastic processes. Intuitively, if two
functions behave very differently, forcing them to share the same basis
functions would affect adversely the approximation accuracy. The second
problem is the balancing issue of different terms in the loss function in MTL.
It is shown in the literature [29] that PINNs, trained in single-task
learning, are already deeply influenced by the weights in front of different
terms in the loss function, e.g., data loss, boundary condition loss, PDE
residual loss. This issue may be more complex in training MH-PINNs, because in
MTL the loss function is commonly defined as weighted summation of task-
specific loss. The last major problem is MH-PINNs for synergistic learning. In
this paper, we only studied one example in function approximation and
presented empirical evidence. More work for the understanding of synergistic
learning with MH-PINNs along both the theoretical and computational directions
should be pursued in the future.

## Acknowledgments

We would like to thank Professor Xuhui Meng of Huazhong University of Science
and Technology for helpful discussions. This work was supported by: the
Vannevar Bush Faculty Fellowship award (GEK) from ONR (N00014-22-1-2795); the
U.S. Department of Energy, Advanced Scientific Computing Research program,
under the Scalable, Efficient and Accelerated Causal Reasoning Operators,
Graphs and Spikes for Earth and Embedded Systems (SEA-CROGS) project, DE-
SC0023191; and by the MURI/AFOSR FA9550-20-1-0358 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. J. Ablowitz and A. Zeppetella, Explicit solutions of Fisher’s equation for a special wave speed, Bulletin of Mathematical Biology, 41 (1979), pp. 835–840.
  * [3] B. Bahmani and W. Sun, Training multi-objective/multi-task collocation physics-informed neural network with student/teachers transfer learnings, arXiv preprint arXiv:2107.11496, (2021).
  * [4] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, Model reduction and neural networks for parametric PDEs, arXiv preprint arXiv:2005.03180, (2020).
  * [5] R. Caruana, Multitask Learning: A Knowledge-Based Source of Inductive Bias, in International Conference on Machine Learning, 1993.
  * [6] X. Chen, J. Duan, and G. E. Karniadakis, Learning and meta-learning of stochastic advection–diffusion–reaction systems from sparse measurements, European Journal of Applied Mathematics, 32 (2021), pp. 397–420.
  * [7] S. Desai, M. Mattheakis, H. Joy, P. Protopapas, and S. Roberts, One-shot transfer learning of physics-informed neural networks, arXiv preprint arXiv:2110.11286, (2021).
  * [8] J. V. Dillon, I. Langmore, D. Tran, E. Brevdo, S. Vasudevan, D. Moore, B. Patton, A. Alemi, M. Hoffman, and R. A. Saurous, TensorFlow Distributions, arXiv preprint arXiv:1711.10604, (2017).
  * [9] L. Dinh, J. Sohl-Dickstein, and S. Bengio, Density Estimation Using Real NVP, arXiv preprint arXiv:1605.08803, (2016).
  * [10] C. Finn, P. Abbeel, and S. Levine, Model-agnostic meta-learning for fast adaptation of deep networks, in International Conference on Machine Learning, PMLR, 2017, pp. 1126–1135.
  * [11] M. Germain, K. Gregor, I. Murray, and H. Larochelle, MADE: Masked autoencoder for distribution estimation, in International Conference on Machine Learning, PMLR, 2015, pp. 881–889.
  * [12] X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, 2010, pp. 249–256.
  * [13] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, Generative adversarial networks, Communications of the ACM, 63 (2020), pp. 139–144.
  * [14] S. Goswami, K. Kontolati, M. D. Shields, and G. E. Karniadakis, Deep transfer operator learning for partial differential equations under conditional shift, Nature Machine Intelligence, 4 (2022), pp. 1155–1164.
  * [15] L. Guo, H. Wu, and T. Zhou, Normalizing field flows: Solving forward and inverse stochastic differential equations using physics-informed flow models, Journal of Computational Physics, 461 (2022), p. 111202.
  * [16] J. Ho, A. Jain, and P. Abbeel, Denoising diffusion probabilistic models, Advances in Neural Information Processing Systems, 33 (2020), pp. 6840–6851.
  * [17] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nature Reviews Physics, 3 (2021), pp. 422–440.
  * [18] R. E. Kass, The validity of posterior expansions based on Laplace’s method, Bayesian and Likelihood Methods in Statistics and Econometrics, (1990), pp. 473–488.
  * [19] Y. Khoo, J. Lu, and L. Ying, Solving parametric PDE problems with artificial neural networks, European Journal of Applied Mathematics, 32 (2021), pp. 421–435.
  * [20] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling, Improved variational inference with inverse autoregressive flow, Advances in Neural Information Processing Systems, 29 (2016).
  * [21] D. P. Kingma and M. Welling, Auto-encoding variational bayes, arXiv preprint arXiv:1312.6114, (2013).
  * [22] I. Kobyzev, S. J. Prince, and M. A. Brubaker, Normalizing flows: An introduction and review of current methods, IEEE Transactions on Pattern Analysis and Machine Intelligence, 43 (2020), pp. 3964–3979.
  * [23] J. Lao, C. Suter, I. Langmore, C. Chimisov, A. Saxena, P. Sountsov, D. Moore, R. A. Saurous, M. D. Hoffman, and J. V. Dillon, tfp. mcmc: Modern Markov chain Monte Carlo tools built for modern hardware, arXiv preprint arXiv:2002.01184, (2020).
  * [24] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, arXiv preprint arXiv:2010.08895, (2020).
  * [25] Z. Lin, Z. Zhao, Z. Zhang, H. Baoxing, and J. Yuan, To Learn Effective Features: Understanding the Task-Specific Adaptation of MAML, 2021, https://openreview.net/forum?id=FPpZrRfz6Ss.
  * [26] K. Linka, A. Schafer, X. Meng, Z. Zou, G. E. Karniadakis, and E. Kuhl, Bayesian Physics-Informed Neural Networks for real-world nonlinear dynamical systems, arXiv preprint arXiv:2205.08304, (2022).
  * [27] X. Liu, X. Zhang, W. Peng, W. Zhou, and W. Yao, A novel meta-learning initialization method for physics-informed neural networks, Neural Computing and Applications, (2022), pp. 1–24.
  * [28] 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.
  * [29] L. McClenny and U. Braga-Neto, Self-adaptive physics-informed neural networks using a soft attention mechanism, arXiv preprint arXiv:2009.04544, (2020).
  * [30] 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.
  * [31] R. M. Neal et al., MCMC using Hamiltonian dynamics, Handbook of Markov Chain Monte Carlo, 2 (2011), p. 2.
  * [32] G. Papamakarios, E. T. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, Normalizing Flows for Probabilistic Modeling and Inference., J. Mach. Learn. Res., 22 (2021), pp. 1–64.
  * [33] G. Papamakarios, T. Pavlakou, and I. Murray, Masked autoregressive flow for density estimation, Advances in Neural Information Processing Systems, 30 (2017).
  * [34] M. Penwarden, S. Zhe, A. Narayan, and R. M. Kirby, Physics-informed neural networks (PINNs) for parameterized PDEs: A metalearning approach, arXiv preprint arXiv:2110.13361, (2021).
  * [35] P. Perdikaris, D. Venturi, and G. E. Karniadakis, Multifidelity information fusion algorithms for high-dimensional systems and massive data sets, SIAM Journal on Scientific Computing, 38 (2016), pp. B521–B538.
  * [36] 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).
  * [37] A. Raghu, M. Raghu, S. Bengio, and O. Vinyals, Rapid Learning or Feature Reuse? Towards Understanding the Effectiveness of MAML, arXiv preprint arXiv:1909.09157, (2019).
  * [38] 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.
  * [39] S. Ruder, An overview of multi-task learning in deep neural networks, arXiv preprint arXiv:1706.05098, (2017).
  * [40] P. Thanasutives, M. Numao, and K.-i. Fukui, Adversarial multi-task learning enhanced physics-informed neural networks for solving partial differential equations, in 2021 International Joint Conference on Neural Networks (IJCNN), IEEE, 2021, pp. 1–9.
  * [41] H. Wang, H. Zhao, and B. Li, Bridging multi-task learning and meta-learning: Towards efficient training and effective adaptation, in International Conference on Machine Learning, PMLR, 2021, pp. 10991–11002.
  * [42] 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.
  * [43] 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.
  * [44] 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.
  * [45] M. Yang and J. T. Foster, Multi-output physics-informed neural networks for forward and inverse PDE problems with uncertainties, Computer Methods in Applied Mechanics and Engineering, (2022), p. 115041.
  * [46] W. Zhong and H. Meidani, PI-VAE: Physics-Informed Variational Auto-Encoder for stochastic differential equations, arXiv preprint arXiv:2203.11363, (2022).
  * [47] Z. Zou, X. Meng, A. F. Psaros, and G. E. Karniadakis, NeuralUQ: A comprehensive library for uncertainty quantification in neural differential equations and operators, arXiv preprint arXiv:2208.11866, (2022).

##  Appendix A Details of NN architectures and training hyperparameters

For all examples in Secs. 4 and 5, MH-PINNs are implemented as fully-connected
NNs (FNNs) with 3 nonlinear hidden layers, each of which is equipped with $50$
neurons and hyperbolic tangent activation function. The number of heads is the
same as the number of tasks in the corresponding examples: $1,000$ in Sec.
4.1, $2,000$ in Secs. 4.2, 4.3 and 4.4, and $10,000$ in Sec. 4.5. Weights in
the body of MH-PINNs are initialized with Glorot uniform initialization [12]
and biases are initialized with zero, while heads are initialized by sampling
from random normal distribution with $0.05$ standard deviation, for fast
training of NFs and better performance of the learned generators.

Except for the forward problem in Sec. 4.2, NFs in this paper are chosen to be
MAF [33] with $10$ bijectors, i.e. the invertible map in NFs, each of which is
a MADE [11], a NN with masked dense layers, with two nonlinear hidden layers
equipped with $100$ neurons and ReLU activation function. The RealNVP [9] and
IAF [20] used in the forward problem in Sec. 4.2 also have $10$ bijectors,
each of which is a NN with two nonlinear hidden layers equipped with $100$
neurons and ReLU activation function. The implementation mostly follows the
instructions of TensorFlow Probability library [8] for NFs.

PI-GANs [44] implemented in Sec. 4.2 have the following architecture: the
discriminator is a FNN with 3 nonlinear hidden layers, each of which is
equipped with $128$ neurons and Leaky ReLU activation function; the generator
that takes as input $t$ is a FNN with 3 nonlinear hidden layers, each of which
is equipped with $50$ neurons and hyperbolic tangent activation function; the
other generator takes as input a Gaussian random variable in $50$ dimensions
with zero mean and identity covariance matrix, and is implemented as a FNN
with 3 nonlinear hidden layers, each of which has $128$ neurons and hyperbolic
tangent activation function. The input dimensions of those 3 FNNs are $65$,
$1$ and $50$, and the output dimensions are $1$, $50$, $50$, respectively.

For the training of MH-PINNs, full-batch training is deployed with Adam
optimizer for $50,000$ iterations. For the training of NFs, except for the
forward problem in Sec. 4.2, mini-batch training is deployed with batch size
being $100$ and Adam optimizer for $1,000$ epochs. NFs in the forward problem
in Sec. 4.2 are trained for $500$ epochs instead, and $L_{2}$ regularization
is imposed to the parameters of RealNVP for better performance. For all NFs,
to achieve stable training, a hyperbolic tangent function is imposed on the
logarithm of the scale, computed from each bijector, such that the logarithm
of the scale lies in $(-1,1)$. For the training of PI-GANs, min-batch training
is deployed with batch size being $100$ and Adam optimizer for $100,000$
iterations. Besides, the same as in [44, 30], physics-informed Wasserstein
GANs (PI-WGANs) with gradient penalty are employed, in which the coefficient
for gradient penalty is set to be $0.1$. Iteratively, 5 updates of the
discriminator are performed and followed by 1 update of the generators. Except
in training PI-GANs, the learning rate of Adam optimizer is set to be
$10^{-3}$ and other hyperparameters of Adam are set as default. In training
PI-GANs, the learning rate is set to be $10^{-4}$, $\beta_{1}=0.5$ and
$\beta_{2}=0.9$ in Adam optimizer for both discriminator and generators.

Training of MH-PINNs, NFs, and PI-GANs was all performed on a single NVIDIA
TITAN Xp GPU. The L-HYDRA code for TensorFlow implementation along with some
representative examples will be released on GitHub once the paper is accepted.

##  Appendix B Details for performing Bayesian inference

Hamiltonian Monte Carlo (HMC) [31] is employed in all Bayesian inference
examples for uncertainty quantification (UQ) while Laplace approximation [18]
is only employed in the first example. In this paper, HMC with adaptive step
size [23] is used, in which the initial step size is set to be either $0.1$ or
$0.01$, tuned for better acceptance rate. The number of burn-in samples and
the number of posterior samples are set to be $1,000$. The number of steps for
the leapfrog scheme is set to be either $30$ or $50$, also tuned for better
acceptance rate. NeuralUQ library [47] for UQ in scientific machine learning
is used as a tool for physics-informed Bayesian inference in the downstream
tasks. The ideal acceptance rate in HMC, as discussed in [30, 47], is around
$60\%$. In this paper, we found chains with acceptance rate from $50\%$ to
$80\%$ acceptable.