# Simulator-free Solution of High-dimensional Stochastic Elliptic Partial
Differential Equations using Deep Neural Networks

Sharmila Karumuri<EMAIL_ADDRESS>Rohit Tripathy<EMAIL_ADDRESS>Ilias Bilionis<EMAIL_ADDRESS>[ Jitesh Panchal<EMAIL_ADDRESS>School
of Mechanical Engineering,
585 Purdue Mall, Purdue University, West Lafayette, IN 47907-2088, USA

###### Abstract

Stochastic partial differential equations (SPDEs) are ubiquitous in
engineering and computational sciences. The stochasticity arises as a
consequence of uncertainty in input parameters, constitutive relations,
initial/boundary conditions, etc. Because of these functional uncertainties,
the stochastic parameter space is often high-dimensional, requiring hundreds,
or even thousands, of parameters to describe it. This poses an insurmountable
challenge to response surface modeling since the number of forward model
evaluations needed to construct an accurate surrogate grows exponentially with
the dimension of the uncertain parameter space; a phenomenon referred to as
the curse of dimensionality. State-of-the-art methods for high-dimensional
uncertainty propagation seek to alleviate the curse of dimensionality by
performing dimensionality reduction in the uncertain parameter space. However,
one still needs to perform forward model evaluations that potentially carry a
very high computational burden. We propose a novel methodology for high-
dimensional uncertainty propagation of elliptic SPDEs which lifts the
requirement for a deterministic forward solver. Our approach is as follows. We
parameterize the solution of the elliptic SPDE using a deep residual network
(ResNet). In a departure from traditional squared residual (SR) based loss
function for training the ResNet, we introduce a physics-informed loss
function derived from variational principles. Specifically, our loss function
is the expectation of the energy functional of the PDE over the stochastic
variables. We demonstrate our solver-free approach through various examples
where the elliptic SPDE is subjected to different types of high-dimensional
input uncertainties. Also, we solve high-dimensional uncertainty propagation
and inverse problems.

###### keywords:

stochastic elliptic partial differential equations , curse of dimensionality ,
deep neural networks , residual networks , energy functional , physics-
informed , high-dimensional uncertainty propagation , inverse problems

††journal: Journal of Computational Physics

url]https://www.predictivesciencelab.org/

##  1 Introduction

With the rapid increase in computing resources [1], numerical methods for the
solution of partial differential equations (PDEs) that govern physical systems
have become an integral part of modern computational science [2]. In many
engineering/scientific applications of interest, inputs to the governing PDEs
are unknown exactly. Such uncertainties are modeled using the framework of
probability theory, thereby giving rise to stochastic partial differential
equations (SPDEs). The uncertainty in input data can arise from multiple
sources - unknown (or partially known) material properties / constitutive
relationships, external loads, initial conditions (ICs), boundary conditions
(BCs), physical geometry (arising from manufacturing imperfections), etc.
Naturally, this leads to the question of how to ensure reliable and robust
predictions of the behavior of the physical system under consideration.
Answering this question is at the heart of research efforts in the field of
uncertainty quantification (UQ) [3, 4]. In particular, the probabilistic
assessment of the effect of input uncertainties on output quantities of
interest (QoI) is known as the forward UQ or uncertainty propagation (UP)
problem.

The most straightforward approach to tackling the UP problem is the Monte
Carlo (MC) method [5]. The MC approach can be summarized as follows - one can
obtain estimates of the statistics of QoIs by computing averages over random
samples. The variance of the MC estimate vanishes in the limit of infinite
samples. A remarkable feature of the MC method is that the statistical
convergence of the MC estimate is independent of the number of the stochastic
dimensions. This makes MC method a highly attractive tool for numerical
integration, especially in high dimensions. Consequently, the MC method and
it’s advanced variants have been extensively applied to a variety of UQ
problems [6, 7, 8, 9, 10]. Unfortunately, the number of samples needed to
obtain convergent statistics from the MC method is typically large. This makes
the application of MC unfeasible for UP in sophisticated modern PDE solvers
because of the high computational cost associated with generating each
individual sample of the QoI.

The standard approach to dealing with expensive numerical PDE solvers is to
replace them with a cheap-to-evaluate surrogate model (or response surface).
The idea of the surrogate approach is that one can utilize information
collected from a finite number of runs of the expensive PDE solver on
carefully selected (and potentially adaptive) input locations, to construct a
map that links the input uncertain parameters to the QoIs. Since the surrogate
is cheap to evaluate, one can apply standard MC methods to propagate
uncertainties through the PDE. This approach to tackling the UP problem has
been applied with great success on a very diverse set of applications with
moderate stochastic dimensionality, such as fracture mechanics [11],
biological systems [12], molecular dynamics [13], nuclear engineering [14],
etc. Traditional choices for surrogate models include Gaussian process
regression (GPR) (or Kriging) [15, 16, 17, 18, 19], polynomial chaos expansion
[20, 21, 22, 23], and radial basis functions [24, 25].

Inspite of the indisputable success of traditional surrogate models on tasks
with low/moderate stochastic dimensionality, approaches such as Gaussian
processes and polynomial chaos have not been successfully scaled to high
dimensions. The task of constructing a response surface becomes exponentially
difficult as the number of input dimensions increases, a phenomenon widely
referred to as the curse of dimensionality [26]. The implication of the curse
of dimensionality is, essentially, that as the input dimensionality grows
linearly, one has to perform an exponentially increasing number of forward
model evaluations to maintain the accuracy of the response surface. In fact,
if one considers the task of approximating a generic nonlinear, multivariate
function with scalar response, the computational time needed for performing
sufficient forward model evaluations quickly becomes unfeasible even for
inexpensive computer codes [27].

Given that the naive construction of response surfaces with high stochastic
dimensionality is a futile approach, a typical workaround is to perform
dimensionality reduction on the stochastic parameter space. The simplest
approach to dimensionality reduction involves a ranking of the input
dimensions in order of importance and then rejecting ‘unimportant’ dimensions.
This is the approach adopted by methods such as sensitivity analysis [28] and
automatic relevance determination [29]. The most common approach to
dimensionality reduction involves projecting stochastic inputs onto a low-
dimensional manifold. In applications characterized by functional
uncertainties (such as flow through porous media), the infinite dimensional
uncertainty is reduced to finite dimensions through the celebrated Karhunen-
Loéve expansion (KLE) [30]. The KLE involves computing the eigendecomposition
of the covariance function of the uncertain parameter and using the decay of
the eigenvalues to approximate the infinite dimensional uncertainty as a
linear combination of the orthogonal eigenfunctions (corresponding to the
retained eigenvalues). In essence this procedure is a linear projection of the
input uncertainty onto a finite dimensional vector subspace. Analogously, when
the uncertain parameters are finite dimensional, a linear projection is
performed on the basis of the eigendecomposition of an empirical covariance
matrix. In the machine learning (ML) community, this is commonly referred to
as the principal component analysis (PCA) [31]. Although KLE and PCA have been
applied successfully to numerous applications, they overestimate the intrinsic
dimensionality because of the fact they (i) recover _only_ linear manifolds in
the input space, and (ii) do not take into consideration information in the
model outputs. Kernel principal component analysis (KPCA) [32, 33] alleviates
the first drawback of the aforementioned techniques by performing the
eigendecomposition on a high (potentially infinite) dimensional space obtained
through a nonlinear transformation of the original inputs. The recently
popularized method of active subspaces (AS) [34, 35, 36, 37, 38, 39, 40], on
the other hand, performs linear dimensionality reduction by performing
eigendecomposition on an empirical covariance matrix constructed from samples
of the model output gradients - thereby alleviating the second drawback. To
bypass the necessity of obtaining gradients (often unfeasible for
sophisticated PDE solvers), [19] proposed a methodology wherein the orthogonal
AS transformation is subsumed into the covariance kernel of GPR. Finally,
recent work from [41] overcomes both limitations, by using DNNs to generalize
the gradient-free approach to AS by recovering nonlinear manifolds in the
input space. At the same time, [42] recast the surrogate modeling task as an
image-to-image regression problem, mapping a snapshot of the input uncertainty
to a snapshot of the PDE solution. By leveraging recent advances in deep
convolutional networks, the authors demonstrate this approach on challenging
high-dimensional surrogate modeling tasks in heterogeneous media [42, 43].

While rapid strides have been made in the developing techniques for high-
dimensional surrogate models, a fundamental limitation of existing
methodologies is that one must still perform repeated evaluations of the
forward PDE solver. Dimensionality reduction in the stochastic parameter space
and adaptive design of experiments can only take one so far. In this work, we
approach the task of UP through SPDEs with high-dimensional uncertainties with
a novel paradigm freed from the shackles of deterministic PDE solvers. The
summary of our approach is that we seek to approximate the field variables in
SPDEs as a parameterized function of all relevant input parameters. These
include the stochastic parameters, spatial/temporal, etc. We seek a flexible
parameterization of the field variables which can accurately approximate
complex nonlinear maps and has closed-form gradients with respect to the input
variables. Naturally, the function approximator of choice to represent the
solution of SPDE is DNNs; specifically we use a deep residual network [44, 45,
46]. DNNs are a class of highly flexible and scalable nonlinear function
approximators [47]. It is well-known that under mild conditions, neural
networks with a single hidden layer are universal function approximators [48].
While the idea of DNNs is not novel, their usage in practical applications is
rather recent - thanks in large part to the widespread availability of cheap
computing resources, the active development and maintainence of powerful
automatic differentiation (AD) capable libraries such as Tensorflow [49],
PyTorch [50], MxNet [51], and theoretical advances in stochastic optimization
[52, 53, 54]. Furthermore, in a departure from squared residual (SR) based
minimization used in related works [55, 56], we leverage the variational
principle associated with the analogous deterministic PDE and formulate an
energy functional for the SPDE by computing an expectation over the stochastic
parameters. The SPDE is then solved by minimizing the energy functional with
respect to the DNN parameters.

The energy functional approach has been previously introduced in other works.
For example, in [57] the authors solved high-dimensional deterministic PDEs
from variational principles using DNN approximators, and in [58] they solved
SPDEs with input random fields up to $100$ dimensions following a solver free
approach. Another example is [59] in which the authors build DNN approximators
for SPDEs with convolutional networks. An exhaustive literature review of the
above three physics-informed neural-network-based papers for solving elliptic
PDEs and SPDEs is shown in Tab. 1. Our paper adds to the existing literature
in the following way. First, we derive mathematically the variational
principle for elliptic SPDE proving the uniqueness of solution in an
appropriate Hilbert space. Second, we derive an optimization objective that
provably converges to a local minimum in the subset of the Hilbert space
spanned by the DNNs. Third, we numerically demonstrate the benefits of using
the energy functional compared to the integrated squared residual. Fourth, we
train the network for various high-dimensional random input fields, including
non-trivial mixtures of random fields, and we demonstrate good performance
even in out-of-sample validation examples. Finally, we assess the accuracy of
the resulting response surface in high-dimensional uncertainty propagation and
inverse problems.

This manuscript is organized as follows. We begin with a discussion of the
variational formulation of the elliptic SPDE problem in Sec. 2.1 and prove
that the solution of the corresponding stochastic boundary value problem
(SBVP) minimizes an energy functional. In Sec. 2.2, we discuss the formulation
of a DNN approximator for the SBVP solution. We propose a construction for the
trial solution in Sec. 2.2.1 which automatically satisfies the essential (i.e.
Dirichlet) boundary conditions. We then move onto the discussion of deep
residual networks (ResNet) in Sec. 2.2.2 \- our function approximator of
choice for this work. In Sec. 2.2.3, we discuss the minimization of the
proposed energy functional loss over the space of ResNets. Sec. 3, is
dedicated to numerical examples. In Sec. 3.1, we discuss the metrics used for
evaluating the performance of our DNN. In Sec. 3.2, we numerically demonstrate
the benefits of using the energy functional compared to the integrated squared
residual with a 1D SBVP example. In Sec. 3.3, we carry out exhaustive
benchmarking of our methodology by solving a 2D SBVP subjected to various
high-dimensional random input fields. In Sec. 3.4, we also test the
generalizability of our DNN approximators to out-of-distribution input data.
In Sec. 3.5, we solve high-dimensional uncertainty propagation (UP) and
inverse problems using the trained DNN approximators. Finally, Sec. 4 is
devoted to the concluding remarks and an outlook on future work.

Table 1: Detailed literature review of physics-informed DNN based papers for solving elliptic PDEs and SPDEs. Author | Summary | Example problems | DNN type |  | Loss
---
function
Optimizer |  | No. of
---
iterations
Batch size
Weinan E et al., 2017 [57] |  They solve high dimensional deterministic PDEs using energy functional (EF) of the PDE via collocation points in each iteration. |  | Poisson problem
---
(in 2d, 10d, 100d).
| ResNet
---
(10 neurons/layer,
6 layers for 10d),
ResNet
(100 neurons/layer,
6 layers for 100d).
EF | ADAM | 50000 |  | 1000 internal
---
points,
100 boundary
points.
| Poisson problem
---
with Neumann BC
(in 5d,10d).
|  |  |  |
| Transfer learning
---
(in 2d).
|  |  |  |
| Eigen value problem
---
(in 1d, 5d, 10d).
|  |  |  |
Meidani et al., 2019 [58] |  They solve high dimensional SPDEs using both squared residual(SR) and energy functional(EF) forms. |  | Transient diffusion
---
problem smooth
random field (in 100d).
| ResNet
---
(256 neurons/layer,
24 layers)
SR | ADAM | $4.6\times 10^{5}$ | 32 samples
| Transient diffusion
---
problem non-smooth
random field (in 50d).
| ResNet
---
(256 neurons/layer,
24 layers)
SR | ADAM | $7.8\times 10^{5}$ | 32 samples
| Steady state heat
---
eqn. (in 50d).
| ResNet
---
(256 neurons/layer,
20 layers)
EF | ADAM | $1.6\times 10^{5}$ | 32 samples
| Steady state heat eqn.
---
(in 30d) in a spatial
domain with hole.
| ResNet
---
(256 neurons/layer,
20 layers)
EF | ADAM | $1.6\times 10^{6}$ | 32 samples
Zhu et al., 2019 [59] |  | Part 1: Compares convolutional
---
decoder network vs fully
connected dense network for
solving PDEs by solving a
deterministic PDE with fixed
input field.
| a) Input field from
---
GRF KLE1024,
b) Non-linear PDE:
non-linear correction
of darcy law with input
field from GRF KLE1024.
| Convolutional
---
decoder net,
FC-NN (8 layers,
512 neurons/layer).
mixed SR | L-BFGS | 500, 2000 iter. |
| Part 2: Solves SPDE, compares
---
physics constrained surrogate
PCS (without labelled data) vs
data driven surrogate DDS
(with labelled data).
| a) Comparison of PCS
---
and DDS with data from
GRF KLE512,
b) Transfer learning check
on out of distribution
input data.
| Convolutional
---
encoder decoder
network
mixed SR | ADAM |  | 300 epochs
---
for PCS,
200 epochs
for DDS.
8 to 32 images
| Part 3: Builds probabilistic
---
surrogate with reverse KL
formulation.
| a) Input field from
---
GRF KLE100, selects
reference density
parameter beta using
reliability diagram.
Glow | mixed SR | ADAM | 400 epochs | 32 images

##  2 Methodology

###  2.1 Variational formulation of stochastic elliptic partial differential
equations

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, where $\Omega$
is the sample space, $\mathcal{F}$ a $\sigma$-algebra of events, and
$\mathbb{P}$ a probability measure. We follow the standard notation where
upper case letters denote random quantities and lower case letters their
values. Let $\mathcal{X}\subset\mathbb{R}^{d}$ be the spatial domain of
interest with dimension $d=1,2$, or $3$. We assume that $\mathcal{X}$ has a
Lipschitz boundary consisting of two disjoint parts $\Gamma_{D}$ and
$\Gamma_{N}$. We will denote points in $\mathcal{X}$ by $x$ and points in
$\Omega$ by $\omega$.

We are interested in approximating the stochastic process (s.p.)
$U:\mathcal{X}\times\Omega\rightarrow\mathbb{R}$ which solves the stochastic
elliptic PDE (SEPDE):

$-\nabla\cdot\left(A(x,\omega)\nabla
U(x,\omega)\right)+C(x,\omega)U(x,\omega)=F(x,\omega),\;\text{in}\;x\in\mathcal{X},$
(1)

almost surely (a.s) with Dirichlet and Neumann boundary conditions given by

$U(x,\omega)=G_{D}(x,\omega)\;\text{on}\;x\in\Gamma_{D},$ (2)

a.s. and

$n^{T}(x)\left(A(x,\omega)\nabla
U(x,\omega)\right)=G_{N}(x,\omega),\;\text{on}\;x\in\Gamma_{N},$ (3)

a.s., respectively, where $n(x)$ is the normal to the boundary $\Gamma_{N}$ at
$x$ on $\Gamma_{N}$, where $C,F,G_{D},G_{N}$ are scalar s.p.’s, and $A$ is
symmetric $d\times d$ matrix s.p.

We derive a variational principle for the above stochastic boundary value
problem (SBVP) in the form of the following theorem which constitutes the main
result of the paper.

######  Theorem 1.

Variational formulation of SEPDEs. Let $C\in
L^{\infty}(\mathcal{X}\times\Omega)$, $F\in L^{2}(\mathcal{X}\times\Omega)$,
$G_{D}\in L^{2}(\Gamma_{D}\times\Omega)$, $G_{N}\in
L^{2}(\Gamma_{N}\times\Omega)$. Assume that $A\in
L^{2}(\mathcal{X}\times\Omega,\mathbb{R}^{d\times d})$ and that it is a.s.
uniformly elliptic, i.e. there exists an $\alpha>0$ such that

$\sum_{i,j=1}^{d}A_{ij}(x,\omega)v_{i}v_{j}\geq\alpha\parallel
v\parallel^{2},\;\text{for
all}\;v\in\mathbb{R}^{d}\;\text{and}\;x\in\mathcal{X}\;\text{a.s.},$ (4)

and that $C\geq\eta>0$ a.s. Then the SBVP defined by Eqs. (1), (2), and (3)
has a unique solution $U^{*}$ in:

$\begin{array}[]{ccc}\mathcal{U}&:=&\Bigg{\\{}U:\mathcal{X}\times\Omega\rightarrow\mathbb{R}\;\text{s.t.}\;\mathbb{E}\left[\int_{\mathcal{X}}\left(\parallel\nabla
U\parallel^{2}+U^{2}\right)dx\right]<\infty,\\\
&&\text{and}\;U_{|\Gamma_{D}}=G_{D}\;\text{a.s.}\Bigg{\\}},\end{array}$ (5)

and this solution minimizes the functional:

$J[U]:=\mathbb{E}\left[\int_{\mathcal{X}}\left\\{\frac{1}{2}\left(A\nabla
U\cdot\nabla
U+CU^{2}\right)-FU\right\\}dx-\int_{\Gamma_{N}}G_{N}Ud\Gamma_{N}\right].$ (6)

Here $\mathbb{E}$ is an expectation over the stochastic parameter $\omega$.

###### Proof.

The proof consists of three parts. We start by proving that the problem of
minimizing $J[U]$ has a unique solution $U^{*}$ in $\mathcal{U}$. Then, we
show that for almost all $\omega$,  this $U^{*}(\cdot,\omega)$ satisfies the
weak form of the SBVP. Finally, we show that any solution of the weak form of
the SBVP is a stationary point of $J[U]$.

##### Part 1: Existence and uniqueness of variational minimum

Consider the following Hilbert space

$\begin{array}[]{ccc}\mathcal{H}&:=&\Bigg{\\{}U:\mathcal{X}\times\Omega\rightarrow\mathbb{R}\;\text{s.t.}\;\mathbb{E}\left[\int_{\mathcal{X}}\left(\parallel\nabla
U\parallel^{2}+U^{2}\right)dx\right]<\infty\Bigg{\\}},\par\end{array}$ (7)

with inner product

$\langle
U,V\rangle_{\mathcal{H}}:=\mathbb{E}\left[\int_{\mathcal{X}}\left(\nabla
U\cdot\nabla V+UV\right)dx\right],$ (8)

and the corresponding norm

$\parallel U\parallel_{\mathcal{H}}:=\langle
U,U\rangle_{\mathcal{H}}^{\frac{1}{2}}.$ (9)

Notice that

$J[U]=\frac{1}{2}B[U,U]-T[U],$

where

$B[U,V]:=\mathbb{E}\left[\int_{\mathcal{X}}\left(A\nabla U\cdot\nabla
V+CUV\right)dx\right],$ (10)

is a bounded bilinear form and

$T[U]=\mathbb{E}\left[\int_{\mathcal{X}}FUdx+\int_{\Gamma_{N}}G_{N}Ud\Gamma_{N}\right],$
(11)

is a bounded linear functional. Finally, define the test space

$\begin{array}[]{ccc}\mathcal{V}&:=&\Bigg{\\{}V:\mathcal{X}\times\Omega\rightarrow\mathbb{R}\;\text{s.t.}\;\mathbb{E}\left[\int_{\mathcal{X}}\left(\parallel\nabla
V\parallel^{2}+V^{2}\right)dx\right]<\infty,\\\
&&\text{and}\;V_{|\Gamma_{D}}=0\;\text{a.s.}\Bigg{\\}},\end{array}$ (12)

a subspace of $\mathcal{H}$. Notice that $\mathcal{U}$ is an affine subspace
of $\mathcal{H}$, i.e. $\mathcal{U}=\tilde{G}_{D}+\mathcal{V}$ where
$\tilde{G}_{D}$ is a suitable extension of $G_{D}$ to the entire
$\mathcal{X}\times\Omega$.

The fist variation of $J[U]$ with respect to a $V$ in $\mathcal{V}$ is

$\begin{array}[]{ccc}\frac{\delta J[U]}{\delta V}&:=&\lim_{\epsilon\rightarrow
0}\frac{J[U+\epsilon V]-J[V]}{\epsilon}\\\ &=&\lim_{\epsilon\rightarrow
0}\frac{\frac{1}{2}B[U+\epsilon V,U+\epsilon V]-T[U+\epsilon
V]-\frac{1}{2}B[U,U]+T[U]}{\epsilon}\\\ &=&\lim_{\epsilon\rightarrow
0}\frac{\frac{1}{2}B[U,U]+\epsilon
B[V,U]+\frac{1}{2}\epsilon^{2}B[V,V]-T[U]-\epsilon
T[V]-\frac{1}{2}B[U,U]+T[U]}{\epsilon}\\\ &=&B[V,U]-T[V].\end{array}$ (13)

Since any local minimum of $J[U]$ must have a vanishing first variation, it
suffices to show that the problem of finding a $U\in\mathcal{U}$ such that

$B[V,U]=T[V],$ (14)

for all $V\in\mathcal{V}$, has a unique solution. This is guaranteed by the
Lax-Milgram theorem if we show that $B$ is coercive. Indeed, we have:

$\begin{array}[]{ccc}B[V,V]&=&\mathbb{E}\left[\int_{\mathcal{X}}\left(A\nabla
V\cdot\nabla V+CV^{2}\right)dx\right]\\\
&\geq&\mathbb{E}\left[\int_{\mathcal{X}}\left(\alpha\parallel\nabla
V\parallel^{2}+\eta V^{2}\right)dx\right]\\\
&\geq&\min\\{\alpha,\eta\\}\parallel V\parallel_{\mathcal{H}}^{2}.\end{array}$
(15)

This concludes the proof that $J[U]$ has a unique minimum
$U^{*}\in\mathcal{U}$.

##### Part 2: The extremum $U^{*}$ solves a.s. the weak form of the BVP

Pick $V(x,\omega)=v(x)\Phi(\omega)$ for a $v$ in the Sobolev space
$H^{1}_{0}(\mathcal{X},\Gamma_{D})$ and $\Phi\in L^{2}(\Omega)$. Eq. (14)
implies that:

$\displaystyle B[v\Phi,U^{*}]$ $\displaystyle=$ $\displaystyle T[v\Phi]$
$\displaystyle\Rightarrow\mathbb{E}\left[\int_{\mathcal{X}}\left(A\nabla
U^{*}\cdot\nabla(v\Phi)+CU^{*}(v\Phi)\right)dx\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\int_{\mathcal{X}}F(v\Phi)dx+\int_{\Gamma_{N}}G_{N}(v\Phi)d\Gamma_{N}\right]$
$\displaystyle\Rightarrow\mathbb{E}\left[\Phi\int_{\mathcal{X}}\left(A\nabla
U^{*}\cdot\nabla v+CU^{*}v\right)dx\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\Phi\left\\{\int_{\mathcal{X}}Fvdx+\int_{\Gamma_{N}}G_{N}vd\Gamma_{N}\right\\}\right],$

and since $\Phi(\omega)$ is arbitrary:

$\int_{\mathcal{X}}\left(A\nabla U^{*}\cdot\nabla
v+CU^{*}v\right)dx=\int_{\mathcal{X}}Fvdx+\int_{\Gamma_{N}}G_{N}vd\Gamma_{N}\;\text{a.s.},$
(16)

i.e. $U^{*}(\cdot,\omega)$ satisfies the weak formulation of the SBVP a.s.

##### Part 3: Any weak solution of the SBVP solves the variational problem

This is readily seen by following the arguments of Part 2 backwards. ∎

######  Remark 1.

Note that it is also possible to prove the theorem for $C=0$. To prove this,
notice that:

$\begin{array}[]{ccc}B[V,V]&=&\mathbb{E}\left[\int_{\mathcal{X}}\left(A\nabla
V\cdot\nabla V\right)dx\right]\\\
&\geq&\mathbb{E}\left[\int_{\mathcal{X}}\left(\alpha\parallel\nabla
V\parallel^{2}\right)dx\right]\\\
&=&\alpha\mathbb{E}\left[\int_{\mathcal{X}}\parallel\nabla
V\parallel^{2}dx\right]\\\
&\geq&\alpha\rho\mathbb{E}\left[\int_{\mathcal{X}}\left(\parallel\nabla
V\parallel^{2}+V^{2}\right)dx\right]\\\ &=&\alpha\rho\parallel
V\parallel_{\mathcal{H}}^{2},\end{array}$ (17)

for some $\rho>0$ that depends only on the Lebesgue measure of $\mathcal{X}$
and on its dimension. Here, going from the third to the fourth step we used
the fact that the $H^{1}$ seminorm and the $H^{1}$ norm are equivalent in
$H^{1}_{0}(\mathcal{X})$, see Chapter 1 of [60].

###  2.2 Building a DNN approximator

From this point on, we assume that the coefficient $C$, source $F$ and the
boundary values $G_{D}$, $G_{N}$ are not stochastic. Also, we assume that the
random field $A(x,\omega)$ is a $d\times d$ scalar matrix with the diagonal
element as $\widetilde{A}(x,\omega)$. Hence Eq. (1), (2) and (3) reduces to
the following,

$\begin{split}-\nabla\cdot\left(\widetilde{A}(x,\omega)\nabla
U(x,\omega)\right)+C(x)U(x,\omega)=F(x),\;\text{in}\;x\in\mathcal{X},\\\\[1.0pt]
U(x,\omega)=G_{D}(x)\;\text{on}\;x\in\Gamma_{D},\\\\[1.0pt]
n^{T}(x)\left(\widetilde{A}(x,\omega)\nabla
U(x,\omega)\right)=G_{N}(x),\;\text{on}\;x\in\Gamma_{N}.\end{split}$ (18)

We assume that the input random field $\widetilde{A}$ can be characterized
with finite number of variables i.e. say
$\widetilde{A}(x,\omega)\approx\widetilde{A}(x,\Xi(\omega))$ where
$\Xi(\omega)$ is the random vector representing the flattened image(or
discretized version) of the input random field $\widetilde{A}$ in
$\mathbb{R}^{d_{\xi}}$ and ${d_{\xi}}$ is number of pixels (or number of
discretizations). We consider $\widetilde{A}(x,\Xi(\omega))$ to be a piecewise
constant function over the pixels of the image. $U(x,\Xi(\omega))$ is the
corresponding solution response i.e. field of interest value corresponding to
the associated spatial location and image of the input random field.

In this work, given samples of the input field,
$\bigg{\\{}\Xi^{(k)}(\omega)\bigg{\\}}^{N}_{k=1}$, ${\Xi^{(k)}(\omega)}\sim
p(\Xi(\omega))$ we wish to learn the solution of the SPDE where
$p(\Xi(\omega))$ is the corresponding probability density assumed or learned
from the data of, say, images of permeability fields, conductivity fields,
micro-structures etc.

####  2.2.1 Extension to boundary conditions:

The solution $U(x,\Xi(\omega))$ is obtained by minimizing the functional $J$
in Eq. (6). Note that Neumann boundary conditions $G_{N}$ are already included
in the functional $J$, but the Dirichlet boundary conditions $G_{D}$ (see Eq.
(5)) have to be imposed as constraints leading to a constrained optimization
problem.

The resulting constrained optimization problem may be addressed in a number of
ways. One way to impose the Dirichlet condition on Eq. (6) is by using
constrained optimization methods such as penalty formulations [55, 57],
Lagrange multipliers, or active set methods [61]. Another way is to model
$U(x,\Xi(\omega))$ in such a way that the constraints are automatically
satisfied by construction, thus changing the original constrained optimization
problem to an unconstrained one [62]. We resort to the latter way by writing
$U(x,\Xi(\omega))$ as the sum of two terms [62],

$\displaystyle U(x,\Xi(\omega))=B(x)+K(x)N(x,\Xi(\omega)),$ (19)

where $N(x,\Xi(\omega))$ is approximated using a deep residual network (see
Sec. 2.2.2) with inputs as $x$ and $\Xi(\omega)$. The term $B(x)$ contains no
adjustable parameters and ensures that the Dirichlet boundary conditions are
satisfied. The second term $K(x)$ is constructed such that it becomes zero on
the Dirichlet boundary. $B$ and $K$ can also be formulated as DNNs given
boundary data [63]. In this work, we limit ourselves to a manual construction
of $B$ and $K$ without using DNNs.

####  2.2.2 The space of deep residual neural networks

We represent $N(x,\Xi(\omega))$ with a DNN. In particular, $N(x,\Xi(\omega))$
is chosen to be a deep residual network or ‘ResNet’ [44]. It has been
empirically demonstrated in numerous works that adding residual blocks is an
effective method for combating the problem of vanishing gradients and
consequently allow robust training of deeper networks [64, 65, 66].
Furthermore, deep ResNets have, recently, been sucessfully applied to the task
of data driven approximation of governing equations of dynamical systems in
[46].

(a)

(b)

Figure 1: 1(a) schematic of our deep ResNet - initial dense layer followed by
$K$ residual blocks and an output linear layer; 1(b) a single residual block
with $L$ layers each having $n$ neurons.

We denote the approximation of $N(x,\Xi(\omega))$ as
$\widehat{N}(x,\Xi(\omega);\theta)$ where ‘$\theta$’ are parameters
(i.e.weights and biases) of the ResNet. The network accepts as input, $x$ and
$\Xi(\omega)$, and produces the output $\widehat{N}(x,\Xi(\omega);\theta)$,
i.e. it defines a map from $\mathbb{R}^{{m}}$, $m=d+d_{\xi}$, to $\mathbb{R}$.

The structure of our deep ResNet is as follows. We begin with a ‘dense
layer’(or fully connected layer) which performs a linear transformation of the
$m$-dimensional input into a $n$-dimensional space. The output of the
computation from this initial dense layer is passed through $K$ ‘residual’
blocks, each with $L$ dense layers having ‘$n$’ neurons each and with
activation function imposed. Finally the output from this computation is
passed through a dense layer which performs linear transformation of the
$n$-dimensional output to a single value $\widehat{N}(x,\Xi(\omega);\theta)$.
Fig. 1(a) shows complete schematic of our deep ResNet.

The residual blocks differentiate classic DNNs from ResNets. A residual block
comprises of an additional residual connection between its input and output,
i.e. the original input to a residual block is added to the block activation
before passing on to the next stage of computation in the network. A schematic
of a residual block is shown in the Fig. 1(b). Of the numerous choices for the
nonlinear activation function in our deep networks, we elect to use the so-
called ‘Swish’ activation function, which is defined as follows:

$\displaystyle\sigma(z)=zS(\beta z)=\frac{z}{1+e^{-(\beta z)}},$ (20)

where, $S$ is the sigmoid activation function and $\beta$ is either a user-
defined or a tunable parameter. The Swish function was introduced in [67], who
empirically demonstrated the superiority of this activation function compared
to standard rectified linear unit or ReLU activation function, for training
very deep networks. In this work, $\beta$ is set to 1.

Mathematically, we start with the inputs $x$ and $\Xi(\omega)$ and we have:

$z^{(1,0)}=W^{(1,0)}_{0}x+W^{(1,0)}_{1}\Xi(\omega)+b^{(1,0)},$ (21)

followed by

$z^{(i,1)}:=\sigma\left(W^{(i,1)}z^{(i-1,0)}+b^{(i,1)}\right),\\\ $ (22)
$z^{(i,j)}:=\sigma\left(W^{(i,j)}z^{(i,j-1)}+b^{(i,j)}\right),\\\ $ (23)
$z^{(i,0)}:=z^{(i,L)}+z^{(i-1,0)},$ (24)

for blocks $i=2,\dots,(K+1)$ and for layers $j=2,\dots,L$. Finally, we end
with

$\widehat{N}(x,\Xi(\omega);\theta)=z^{(K+2)}:=w^{(K+2)}z^{(K+1,0)}+b^{(K+2)}.$
(25)

The last term in Eq. (24) is the residual connection which helps to ease the
training of the network [44] and the quantities
$W^{(1,0)}_{0}\in\mathbb{R}^{n\times d}$, $W^{(1,0)}_{1}\in\mathbb{R}^{n\times
d_{\xi}}$, {$W^{(i,1)},W^{(i,j)}\\}\in\mathbb{R}^{n\times n}$,
$w^{(K+2)}\in\mathbb{R}^{1\times n}$ and {$b^{(1,0)},b^{(i,1)},b^{(i,j)}$}
$\in\mathbb{R}^{n\times 1}$, $b^{(K+2)}\in\mathbb{R}$ are the weights and bias
parameters, respectively. Collectively all these are the parameters $\theta$
of our ResNet.

####  2.2.3 Training the ResNet

We plug the approximation $\widehat{N}(x,\Xi(\omega);\theta)$ into Eq. (19)
leading to the following expression for the solution of the SBVP:

$\displaystyle\widehat{U}(x,\Xi(\omega);\theta)=B(x)+K(x)\widehat{N}(x,\Xi(\omega);\theta).$
(26)

Consequently, the energy functional, $J$, is now a function of the free
parameters $\theta$:

$J(\theta):=J[\widehat{U}(x,\Xi(\omega);\theta)].$ (27)

Therefore, the task of training the ResNet is equivalent to solving the
following unconstrained optimization task:

$\theta^{*}=\underset{\theta}{\arg\min}J(\theta).$ (28)

Note that the variational problem obtained in the end i.e. $J(\theta)$ is not
convex despite the fact that the initial problem $J[U]$ is convex.

To proceed, consider the sampling average approximation :

$\displaystyle\hat{J}(\theta):=\frac{|X|}{N_{\xi}}\sum_{k=1}^{N_{\xi}}\Bigg{\\{}\frac{1}{N_{x}}\sum_{i=1}^{N_{x}}\Bigg{[}\frac{1}{2}\Big{[}\widetilde{A}\left(X^{(ki)},\Xi^{(k)}(\omega)\right)\nabla_{x}\widehat{U}\left(X^{(ki)},\Xi^{(k)}(\omega);\theta\right)\cdot\nabla_{x}\widehat{U}\left(X^{(ki)},\Xi^{(k)}(\omega);\theta\right)$
(29)
$\displaystyle+C\left(X^{(ki)}\right)\widehat{U}^{2}\left(X^{(ki)},\Xi^{(k)}(\omega);\theta\right)\Big{]}-F\left(X^{(ki)}\right)\widehat{U}(X^{(ki)},\Xi^{(k)}(\omega);\theta)\Big{)}\Bigg{]}\Bigg{\\}}$
$\displaystyle-\frac{|X_{b}|}{N_{\xi}}\sum_{k=1}^{N_{\xi}}\Bigg{\\{}\frac{1}{N_{b}}\sum_{r=1}^{N_{b}}\Bigg{[}G_{n}\left(X^{(kr)}_{b}\right)\widehat{U}\left(X^{(kr)}_{b},\Xi^{(k)}(\omega);\theta\right)\
\Bigg{]}\Bigg{\\}},$

where $\Xi^{(k)}(\omega),k=1,\dots,N_{\xi}$ are independent identically
distributed (iid) replicas of $\Xi(\omega)$,
$X^{(ki)},i=1,\dots,N_{x},k=1,\dots,N_{\xi}$ are iid random vector’s (r.v.’s)
uniformly distributed in $\mathcal{X}$,
$X^{(kr)}_{b},r=1,\dots,N_{b},k=1,\dots,N_{\xi}$ are iid r.v.’s uniformly
distributed on the Neumann boundary $\Gamma_{n}$, $|X|$ is the Lebesgue
measure of the spatial domain $\mathcal{X}$ and $|X_{b}|$ is the Lebesgue
measure of the Neumann boundary $\Gamma_{n}$. It is trivial to see that
$\hat{J}$ is an unbiased Monte Carlo estimate of the energy functional $J$.

Thus, the ResNet training is recast into a stochastic minimization problem:

$\theta^{*}=\underset{\theta}{\arg\min}\ \mathbb{E}[\hat{J}(\theta)].$ (30)

The unconstrained optimization problem in Eq. (30) is solved through the
adaptive moments (ADAM) optimization method [52], a robust variant of
stochastic gradient descent (SGD) [68]. The ADAM’s update scheme is given by:

$\begin{split}\mathbf{m}_{j+1}&\leftarrow\beta_{1}\mathbf{m}_{j}+(1-\beta_{1})\nabla_{\theta}\hat{J}(\theta_{j}),\\\
\mathbf{v}_{j+1}&\leftarrow\beta_{2}\mathbf{v}_{j}+(1-\beta_{2})[\nabla_{\theta}\hat{J}(\theta_{j})]^{2},\\\
\widetilde{\mathbf{m}}_{j+1}&\leftarrow\frac{\mathbf{m}_{j+1}}{1-\beta_{1}},\\\
\widetilde{\mathbf{v}}_{j+1}&\leftarrow\frac{\mathbf{v}_{j+1}}{1-\beta_{2}},\\\
\theta_{j+1}&\leftarrow\theta_{j}-\alpha\frac{\widetilde{\mathbf{m}}_{j+1}}{\sqrt{\widetilde{\mathbf{v}}_{j+1}}+\epsilon},\end{split}$
(31)

where $\alpha$ is a positive learning rate, $\epsilon$ is a small positive
number used to prevent zero in the denominator, $\beta_{1}$ and $\beta_{2}$
are averaging parameters which are, in principle, tunable. In practice,
default values of $\beta_{1}=0.9$, $\beta_{2}=0.999$, as suggested by [52]
work well and we do not change these quantities.  In the loss function
$\hat{J}(\theta)$ of Eq. (29) we obtain the required derivatives using
automatic differentiation (AD) [45] in TensorFlow [69]. The exact gradient of
the loss function needed for the ADAM update in Eq. (31) is obtained using
backpropagation [70].

##  3 Numerical examples

We begin by discussing metrics for evaluating the performance of our DNN and
we then proceed to an exhaustive comparison study between using an energy
functional loss function vs using an integrated squared residual loss function
for training our ResNets that solve an elliptic SBVP problem. Next, we
demonstrate the performance of our DNN approximator in solving a 2D SBVP
problem for different types of input fields. Then, we show the predictive
ability of our trained DNNs on out-of-distribution inputs. Finally, we solve
UP and inverse problems using our trained DNNs. Our DNNs are implemented using
Keras [71], an application programming interface (API) running on top of
TensorFlow [72]. The code and data used for training and testing will be made
available at https://github.com/PredictiveScienceLab/variational-elliptic-SPDE
upon publication.

###  3.1 Evaluation metrics

The solution predicted by our deep ResNet is compared to the solution obtained
from a finite volume method (FVM) solver implemented in FiPy [73]. To do this,
we discretize the spatial domain into $N_{\text{cells}}$ cells. Input random
fields discretized at the cell centers are fed as input to the FiPy solver.
Then, the solver estimates the corresponding numerical solutions of the SBVP
at these cell locations. We calculate the relative root mean square error
‘$\mathcal{E}$’ between the predicted solution from our DNN and the solution
obtained from a FVM solver for all the $N_{\text{sam}}$ test samples as
follows:

$\mathcal{E}=\sqrt{\frac{\sum_{i=1}^{N_{\text{sam}}}\sum_{j=1}^{N_{\text{cells}}}(U_{i,j}^{\text{FVM}}-U_{i,j}^{\text{DNN}})^{2}}{\sum_{i=1}^{N_{\text{sam}}}\sum_{j=1}^{N_{\text{cells}}}(U_{i,j}^{\text{FVM}})^{2}}},$
(32)

where $U_{i,j}^{\text{FVM}}$ is the FVM solution at the $j^{th}$ cell center
corresponding to the $i^{th}$ sample of the input field, and
$U_{i,j}^{\text{DNN}}$ is the predicted solution from our DNN corresponding to
the same realization of the input field at the same cell center location.

We also evaluate the quality of our deep ResNet’s prediction based on two
metrics:

  1. 1.

The relative $L_{2}$ error metric defined as:

$L_{2}(\mathbf{U}_{\text{DNN}},\mathbf{U}_{\text{FVM}})=\frac{\parallel\mathbf{U}_{\text{FVM}}-\mathbf{U}_{\text{DNN}}\parallel_{2}}{\parallel\mathbf{U}_{\text{FVM}}\parallel_{2}},$
(33)

where, $\|\cdot\|_{\mathrm{2}}$ is the standard Euclidean norm.
$\mathbf{U}_{\text{FVM}}$ and $\mathbf{U}_{\text{DNN}}$ are the FVM solution
vector and the DNN prediction vector corresponding to a particular realization
of the input field.

  2. 2.

The coefficient of determination, (also known as the $R^{2}$ score), defined
as:

$R^{2}=1-\frac{\sum_{k=1}^{N_{cells}}(U_{\text{FVM},k}-U_{\text{DNN},k})^{2}}{\sum_{k=1}^{N_{cells}}(U_{\text{FVM},k}-\bar{U}_{\text{FVM}})^{2}},$
(34)

where, $k$ indexes all the FVM cell centers, $U_{\text{FVM},k}$ and
$U_{\text{DNN},k}$ are the FVM solution and the DNN predicted solution at the
$k^{th}$ cell center respectively, and $\bar{U}_{\text{FVM}}$ is the mean of
$U_{\text{FVM},k}$.

###  3.2 Comparison of Energy Functional (EF) loss function vs integrated
Squared Residual (SR) loss function with a 1D SBVP

We explore the relative merit of using EF based loss function over the
integrated SR based loss function, in training ResNet approximators for
elliptic SBVP solution. For this study, consider an elliptic SBVP in 1D on a
unit length domain:

$-\nabla\cdot(\widetilde{A}(x,\omega)\nabla U(x,\omega))+C\ U(x,\omega)=F,\
\forall x\in\mathcal{X}=[0,1]\subset\mathbb{R}^{1},$ (35)

with $C=15$, $F=10$ and subjected to Dirichlet boundary conditions on both the
ends:

$\begin{split}U=1,\ \forall\ x=0,\\\ U=0,\ \forall\ x=1.\end{split}$ (36)

The stochasticity in Eq. (35) arises from the uncertainty in the spatially-
varying input field $\widetilde{A}(x,\omega)$. The input random field,
$\widetilde{A}$ is modeled as log-normal field -

$\log(\widetilde{A}(x,\omega))\sim\mathrm{GP}(\mu(x),k(x,x^{\prime})),$ (37)

with mean $\mu(x)=0$ and exponential covariance function

$k(x,x^{\prime})=\sigma^{2}exp\left\\{-\frac{|x-x^{\prime}|}{\ell_{x}}\right\\},$
(38)

where, $\ell_{x}$ and $\sigma^{2}$ represent the correlation length/length-
scale and variance of the log input field respectively. We bound this input
field $\widetilde{A}$ uniformly from below by $0.005$ to satisfy stochastic
ellipticity condition and we also bound it uniformly from above by $33$ to
avoid very large input field values.

For this comparison study, we generated $10,000$ bounded train samples and
$1,000$ bounded test samples of this input field
$(\widetilde{A}(x,\omega)\approx(\widetilde{A}(x,\Xi(\omega))$ independently
from length-scale $\ell_{x}$ in $\\{2.0,1.0,0.5,0.1,0.03\\}$ and
$\sigma^{2}=1$ over a uniform grid of $d_{\xi}=100$ cells/pixels.

Following Eq. (19), we account for the Dirichlet boundary conditions in a
hard-way by writing:

$\displaystyle U(x,\Xi(\omega))=(1-x)+(x(1-x))N(x,\Xi(\omega)).$ (39)

We approximate $N$ as $\widehat{N}(x,\Xi(\omega);\theta)$ with our deep ResNet
shown in Fig. 1(a) which accepts $m=d+d_{\xi}=1+100=101$ inputs.

We fix the ResNet architecture and hyper-parameters (see Tab. 2) for each
length-scale and train it with both EF loss (Eq. (29)) and integrated SR loss
which is based on the residual of SPDE in Eq. (35) and is given by (the
notations have same meaning as in Eq. (29))

$\displaystyle\hat{J}_{SR}(\theta):=\frac{|X|}{N_{\xi}}\sum_{k=1}^{N_{\xi}}\Bigg{\\{}\frac{1}{N_{x}}\sum_{i=1}^{N_{x}}\Bigg{[}\Big{[}-\nabla_{x}\cdot\widetilde{A}\left(X^{(ki)},\Xi^{(k)}(\omega)\right)\nabla_{x}\widehat{U}\left(X^{(ki)},\Xi^{(k)}(\omega);\theta\right)$
(40)
$\displaystyle+C\left(X^{(ki)}\right)\widehat{U}\left(X^{(ki)},\Xi^{(k)}(\omega);\theta\right)\Big{]}-F\left(X^{(ki)}\right)\Bigg{]}^{2}\Bigg{\\}}.$

We perform optimization of these loss functions using the ADAM optimizer with
a constant learning rate (see Tab. 2) by randomly picking $N_{\xi}=100$
realizations of the input field from the train samples, and then on each of
picked realizations we uniformly sample $N_{x}=15$ spatial points in
$\mathcal{X}$ for each iteration (here $|X|=1$). Then we estimate the
$\mathcal{E}$ of all the test samples using DNNs trained on both the loss
functions. Training each of these DNN’s took around 25 minutes to 3.4 hours of
computational time depending on the size of the network.

$\ell_{x}$ | $K$ | $L$ | $n$ | $learning\ rate$ |  Number of iterations |  Number of trainable parameters $\theta$
---|---|---|---|---|---|---
2 | 3 | 2 | 10 | 0.001 | 65,000 | 1,691
1 | 3 | 2 | 80 | 0.001 | 65,000 | 47,121
0.5 | 3 | 2 | 130 | 0.001 | 65,000 | 115,571
0.1 | 3 | 2 | 200 | 0.0001 | 70,000 | 261,801
0.03 | 3 | 2 | 400 | 0.0001 | 80,000 | 1,003,601
Table 2: (1D SBVP) ResNet architecture and parameters used for each length-
scale in the comparison study.

We replicate the analysis 6 times with different seeds for each length-scale.
The results are provided in Fig. 2 using box plots. Each box is represented by
the first quartile (bottom line), the median (middle line) and the third
quartile (upper line). The whiskers indicate the variability of the
performance outside the first and third quartiles. The ends of the whiskers
lie at a distance of 1.5 interquartile range from the first/third quartile. We
can clearly see from the figure that both the median and interquartile ranges
of $\mathcal{E}$ are significantly smaller for the DNNs trained with EF loss
over integrated SR loss at each length-scale. We observed that as length-scale
decreases the integrated SR loss fails to train and it just learns the mean-
profile. This study demonstrates the superior behavior of EF loss over
integrated SR loss for training ResNet based approximators of elliptic SBVP.

Figure 2: (1D SBVP) Box plots of relative root mean square error $\mathcal{E}$
based on 6 runs for each length-scales using EF and integrated SR losses.

For the sake of completeness, in Fig. 3 we show the test predictions from the
DNN trained with EF loss function corresponding to input field realizations of
length-scale $0.03$ and variance $1$. We observe that the relative $L_{2}$
error as reported on the headers is less than $0.05$ and the $R^{2}$ score
close to $0.99$, which implies that the predicted solution from DNN matches
truth from FVM very closely. Fig. 4 shows the histograms of the relative
$L_{2}$ errors and $R^{2}$ scores of all the test samples.

Figure 3: (1D SBVP - GRF $\ell_{x}$ $0.03$) Each row corresponds to a randomly
chosen realization of log-input field (left column) from the GRF of length-
scale $\ell_{x}$ $0.03$ test dataset and the corresponding solution response
(right column). The DNN prediction is the red solid line and the FVM solution
is the black dotted line.

(a)

(b)

Figure 4: (1D SBVP - GRF $\ell_{x}$ $0.03$) 4(a) and 4(b) corresponds to
histograms of relative $L_{2}$ errors and $R^{2}$ scores for all the input
field samples in the test dataset of GRF of length-scale $\ell_{x}$ $0.03$.

###  3.3 Stochastic boundary value problem in 2D

Consider the following elliptic SBVP in 2D on a unit square domain -

$-\nabla\cdot(\widetilde{A}(x,\omega)\nabla U(x,\omega))=0,\ \ \forall
x\in\mathcal{X}=[0,1]^{2}\subset\mathbb{R}^{2},$ (41)

with boundary conditions:

$\begin{split}U=1,\ \;\forall\;x_{1}=0,\\\ U=0,\ \;\forall\;x_{1}=1,\\\
n^{T}(x)\left(\widetilde{A}(x,\omega)\nabla U(x,\omega)\right)=0,\
\;\forall\;x_{2}=0\;\text{and}\;x_{2}=1.\end{split}$ (42)

Eq. (41) models steady-state diffusion processes in 2D. The quantity
$\widetilde{A}(x,\omega)$ is a spatially varying diffusion coefficient. The
physical significance of the equation and all terms in it varies from context
to context. For instance, Eq. (41) could be a model for single-phase
groundwater flow, where $\widetilde{A}$ represents the permeability
coefficient and the solution variable $U$ the pressure. Similarly, Eq. (41)
could also be a model for steady-state heat conduction, where $\widetilde{A}$
represents the conductivity coefficient and the solution variable $U$ the
temperature.

We solve this 2D elliptic SPDE in the context of 4 different input fields
($\widetilde{A}$) namely Gaussian Random Field (GRF), warped GRF, channelized
field, and multiple length-scales GRF. Similar to the previous example, all
the input field samples discussed below are constrained uniformly by a lower
bound 0.005 and an upper bound 33.

The first input field dataset is exponential of a GRF, i.e.

$\log(\widetilde{A}(x,\omega))\sim\mathrm{GP}(\mu(x),k(x,x^{\prime})),$ (43)

with the mean function $\mu(x)=0$ and an exponential covariance function (or
kernel)

$k(x,x^{\prime})=\sigma^{2}exp\left\\{-\large\sum_{i=1}^{2}\frac{|x_{i}-x^{\prime}_{i}|}{\ell_{x_{i}}}\right\\}.$
(44)

We set the length-scales and the variance of the field to be $l_{x_{1}}=0.05$,
$l_{x_{2}}=0.08$ and $\sigma^{2}=0.75$ respectively and generate 10,000 train
samples and 2,000 test samples of this input field
$(\widetilde{A}(x,\omega)\approx(\widetilde{A}(x,\Xi(\omega))$ over a uniform
grid of $32\times 32$ cells/pixels.

The second input field dataset considered is an exponential of a warped GRF
(or two-layer GRF), where there are two Gaussian fields and the output of the
first GRF is input to the second GRF. This input field is defined implicitly
through a warped GRF as follows:

$S(x,\omega)\sim GP(\mu(x),k_{1}(x,x^{\prime})),$ (45)

where $S=[S_{1},S_{2}]$,

$\log\widetilde{A}(S,\omega)\sim GP(0,k_{2}(S,S^{\prime})),$ (46)

the mean and covariance functions are chosen as follows:

$\begin{split}\mu(x)=x,\\\\[1.0pt]
k_{1}(x,x^{\prime})=\sigma_{1}^{2}exp\left\\{-\large\sum_{i=1}^{2}\frac{||x_{i}-x^{\prime}_{i}||^{2}}{2\ell_{x}^{2}}\right\\},\\\\[1.0pt]
k_{2}(S,S^{\prime})=\sigma_{2}^{2}exp\left\\{-\large\sum_{i=1}^{2}\frac{||S_{i}-S^{\prime}_{i}||^{2}}{2\ell_{S}^{2}}\right\\},\\\\[1.0pt]
\end{split}$ (47)

with $l_{x}=2$, $l_{S}=0.1$, $\sigma_{1}^{2}=0.25$ and $\sigma_{2}^{2}=0.75$.
We generate 10,000 train samples and 1,000 test samples of this input field
over a uniform grid of $32\times 32$ cells/pixels.

The third input field dataset considered is channelized field [74], defined
with binary values $0.01$ and $1.0$. This dataset is taken from this link
https://github.com/cics-nd/pde-surrogate and these samples are obtained by
cropping $32\times 32$ patches from $1250\times 1250$ image. We have $4,096$
train samples and $512$ test samples of this field.

The fourth type of input field dataset considered is also exponential of a GRF
with exponential kernel and variance $0.75$ similar to the first dataset but
here we lift the assumption of fixed length-scales of the GRF. Instead, we
want our dataset to be made up of realizations of different length-scales. So,
intuitively we would like to have more samples from lower length-scales as
variability in the SPDE solution $U$ will be high for these lower length-
scales. Hence, following the procedure in Algorithm 3 of [41] we obtain $60$
different length-scale pairs (shown in Fig. 5) and for each length-scale pair
we obtain $500$ samples of the input field for training and $150$ samples of
the input field for testing which totals down in all to $60\times 500=30,000$
train samples and $60\times 150=9,000$ test samples over a uniform grid of
$32\times 32$ cells/pixels. We call this dataset as multiple length-scales GRF
dataset from here on. Note that the length-scales shown in Fig. 5 are
constrained by a lower bound of FVM cell size i.e. $(h=\frac{1}{32})$.
Representative samples of these input fields is shown in Fig. 6.

Figure 5: (2D SBVP - Multiple length-scales GRF) Visual representation of
sampled length-scale pairs. Each ‘x’ corresponds to a particular length-scale
sampled.

Figure 6: (2D SBVP) Samples from the 4 different input field datasets i.e.
GRF, warped GRF, channelized field and multiple length-scales GRF over a
uniform grid of $32\times 32$. Log input field samples are shown except the
channelized field.

DNN approximators:  Now we build $4$ DNN approximators using each of the input
fields discussed above. Following Eq. (19), we account for the Dirichlet
boundary conditions in a hard-way by writing:

$\displaystyle U(x,\Xi(\omega))=(1-x_{1})+(x_{1}(1-x_{1}))N(x,\Xi(\omega)).$
(48)

We approximate $N$ as $\widehat{N}(x,\Xi(\omega);\theta)$ with our deep ResNet
shown in Fig. 1(a) which accepts $m=d+d_{\xi}=2+1,024=1,026$ inputs. We
utilize ADAM optimizer with constant learning rate of $0.001$ to optimize the
loss function $\widehat{J}(\theta)$ (Eq. 29) by randomly picking $N_{\xi}=100$
realizations of the input field from the train samples, and then on each of
the picked realizations we uniformly sample $N_{x}=20$ spatial points in
$\mathcal{X}$ for each iteration. Also, $|X|=1$ and $|X_{b}|=1$ here for this
2D example. Tab. 3 shows the parameters of the ResNet architecture for each of
the DNN approximators and the estimates of $\mathcal{E}$ over all the samples
in their respective test datasets. We see that all the DNNs predict well very
with $\mathcal{E}$ less than $5.31\%$. Training these DNN’s took around 5 to 7
hours of computational time depending on the size of the network.

$Datasets$ | $K$ | $L$ | $n$ |  Number of test samples | $\mathcal{E}$ |  Number of trainable parameters $\theta$
---|---|---|---|---|---|---
GRF $\ell_{x}$ $[0.05,0.08]$ | 3 | 2 | 350 | 2,000 | $4.45\%$ | 1,096,901
Warped GRF | 5 | 2 | 300 | 1,000 | $4.68\%$ | 1,211,401
Channelized field | 3 | 2 | 300 | 512 | $5.30\%$ | 850,201
Multiple length-scales GRF | 3 | 2 | 500 | 9,000 | $3.86\%$ | 2,017,001
Table 3: (2D SBVP) ResNet architecture of the 4 DNN approximators and relative
root mean square error $\mathcal{E}$ for test samples with these networks.

Figs. 7 \- 10 show a comparison of the SPDE solution predicted by the 4 DNN
approximators vs solution obtained from the FVM solver for 4 randomly chosen
test samples of their corresponding input fields. We see from the above
figures that DNNs are able to capture fine-scale features in the solution
responses. Fig. 11 shows histograms of relative $L_{2}$ errors and $R^{2}$
scores of the 4 DNN approximators over all the samples in their respective
test dataset.

Figure 7: (2D SBVP - GRF $\ell_{x}$ $[0.05,0.08]$) Each row corresponds to a
randomly chosen realization of log-input field (left column) from the GRF of
length-scales $[0.05,0.08]$ test dataset and the corresponding solution
response from FVM and DNN (middle and right columns). Figure 8: (2D SBVP -
Warped GRF) Each row corresponds to a randomly chosen realization of log-input
field (left column) from the warped GRF test dataset and the corresponding
solution response from FVM and DNN (middle and right columns). Figure 9: (2D
SBVP - Channelized field) Each row corresponds to a randomly chosen
realization of input field (left column) from the channelized field test
dataset and the corresponding solution response from FVM and DNN (middle and
right columns). Figure 10: (2D SBVP - Multiple length-scales GRF) Each row
corresponds to a randomly chosen realization of log-input field (left column)
from the multiple length-scales GRF test dataset and the corresponding
solution response from FVM and DNN (middle and right columns).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 11: (2D SBVP - 4 fields) (11(a)-11(b)), (11(c)-11(d)), (11(e)-11(f))
and (11(g)-11(h)) corresponds to histograms of relative $L_{2}$ errors and
$R^{2}$ scores for all the input field samples in the test datasets of GRF
$\ell_{x}$ $[0.05,0.08]$, warped GRF, channelized field and multiple length-
scales GRF respectively.

We also trained a DNN over all the train samples of the $4$ input field
datasets discussed before. This DNN was trained with same optimizer settings
and batch sizes as the previous $4$ DNN approximators, and has a ResNet
architecture (see Fig. 1(a)) of $K=2$ residual blocks each with $L=2$ layers
having $n=500$ neurons each. The $\mathcal{E}$ over all the test samples in
$4$ datasets comes out to be $4.56\%.$ From Fig. 12 we see that this DNN
(single DNN) is able to predict the solution response for all the input test
datasets very well except the channelized fields. This is likely caused by the
small number of training samples in the channelized field dataset compared to
others. So if an equal proportion of train samples from all the datasets are
given during training then the network performance could be improved further.
Fig. 13 shows comparison plots of DNN and FVM predictions for a few randomly
chosen test samples. Fig. 14 shows histograms of relative $L_{2}$ errors and
$R^{2}$ scores for all the test samples of $4$ input field datasets. From this
study we see that given sufficient amount of train data of all types of fields
a single DNN has the potential to learn all the solution responses.

Figure 12: (2D SBVP - Single DNN) Relative root mean square error,
$\mathcal{E}$ comparison for the 4 input field test datasets between a
separate DNN approximator vs a single DNN approximator.

(a)

(b)

(c)

(d)

Figure 13: (2D SBVP - Single DNN) 13(a), 13(b), 13(c) and 13(d) left-columns
corresponds to realizations of log-input field from the $4$ test datasets -
GRF $\ell_{x}$ $[0.05,0.08]$, warped GRF, channelized field and multiple
length-scales GRF respectively, middle and right columns correspond to
solution response from FVM and DNN trained on all 4 input field train images.

(a)

(b)

Figure 14: (2D SBVP - Single DNN) 14(b) and 14(b) corresponds to histograms of
relative $L_{2}$ errors and $R^{2}$ scores for all the input field samples in
the $4$ test datasets (GRF $\ell_{x}$ $[0.05,0.08]$, warped GRF, channelized
field and multiple length-scales GRF).

###  3.4 Testing DNNs generalizability to out-of-distribution inputs

The generalization capability of our trained DNN approximators in the before
section to out-of-distribution input data/new input distribution data not seen
in training is examined here by considering the following $3$ cases-

Note that similar to the previous example, all the out-of-distribution input
field samples discussed below are also constrained uniformly by a lower bound
0.005 and an upper bound 33.

Case A: Here we study the generalizabilty of DNN approximator trained with GRF
of length-scales $[0.05,0.08]$ dataset (in Sec. 3.3) to 4 sets of out-of-
distribution input data which are GRF with a larger length-scale of
$[0.3,0.4]$ and variance of $0.75$, GRF with a smaller length-scale of
$[0.03,0.04]$ and variance of $0.75$, warped GRF and channelized field. The
relative root mean square error, $\mathcal{E}$ of predicted fields of these 4
out-of-distribution input datasets in shown in Fig. 15. This figure shows that
DNN generalizes very well to all the out-of-distribution inputs except the
channelized field one as expected since it is far from the input field
distribution on which the DNN is trained. Fig. 16 shows few out-of-
distribution input samples and the predicted solution responses from DNN.

Figure 15: (Generalizability test - Case A) Relative root mean square error,
$\mathcal{E}$ of four out-of-distribution input datasets, each with $512$
samples tested on DNN trained with input field samples from GRF of length-
scales $[0.05,0.08]$.

(a)

(b)

(c)

(d)

Figure 16: (Generalizability test - Case A) 16(a), 16(b), 16(c) and 16(d)
left-columns corresponds to realizations of input field from out-of-
distribution input data GRF $\ell_{x}$[0.3,0.4], GRF $\ell_{x}$[0.03,0.04],
warped GRF and channelized field respectively, middle and right columns
correspond to solution response from FVM and DNN trained on GRF data of
length-scales $\ell_{x}$ $[0.05,0.08]$.

Case B: Here we test the generalizabilty of DNN approximator trained with
multiple length-scales GRF dataset (in Sec. 3.3) on new set of realizations
sampled from length-scales not used in training the network. A $10\times 10$
uniform grid of length-scales is generated in the domain $[h,1]^{2}$, and for
each length-scale, 100 samples of the input field are generated i.e. arbitrary
length-scales dataset. The mean of the relative $L_{2}$ errors and mean of the
$R^{2}$ scores for each length-scale pair in this uniform grid is computed and
shown in Fig. 17. Note that, as expected, the accuracy of the DNN decreases
for very fine length-scale values because of large variations in the input
field. Also, Fig. 18 shows DNN predictions for few randomly chosen samples in
the arbitrary length-scales dataset and Fig. 19 shows the histogram of
relative $L_{2}$ errors, $R^{2}$ scores for all the samples in this arbitrary
length-scales dataset. From this study we observe very good generalizability
of our Resnet to predict solutions of SPDE corresponding to even length-scales
not used in its training.

(a)

(b)

Figure 17: (Generalizability test - Case B) 17(a) and 17(b) corresponds to the
mean relative $L_{2}$ errors and mean $R^{2}$ scores of the predicted
solutions for input field samples with arbitrary pairs of length-scales not
used in the DNN training. The black ‘x’ markers correspond to length-scales
used in training the DNN and the solid black dots correspond to arbitrary
length-scales used to test the DNN.

Figure 18: (Generalizability test - Case B) Each row corresponds to a randomly
chosen realization of log-input field (left column) from arbitrary length-
scales data and the corresponding solution response from FVM and DNN trained
on multiple length-scales GRF dataset(middle and right columns).

(a)

(b)

Figure 19: (Generalizability test - Case B) 19(a) and 19(b) corresponds to
histograms of relative $L_{2}$ errors and $R^{2}$ scores for all the input
field samples in the arbitrary length-scales GRF dataset.

Case C: We also test the generalizabilty of DNN approximator trained with
multiple length-scales GRF dataset (in Sec. 3.3) on stratified (or 2-layered)
input fields. Stratified fields are generated by dividing the spatial region
$\mathcal{X}$ into two parts - $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$ through
a line joining the two points $(0,a)$ and $(1,b)$. $a$ and $b$ are sampled
independently from $\mathcal{U}([0,1])$. The stratified input fields are
modeled as follows:

$\log(\widetilde{A}(x,\omega))\sim\sigma^{2}\Big{(}\mathrm{GP}(0,k_{1}(x,x^{\prime}))\mathbb{I}_{\mathcal{X}_{1}}(x)+\mathrm{GP}(0,k_{2}(x,x^{\prime}))\mathbb{I}_{\mathcal{X}_{2}}(x)\Big{)},$
(49)

where, $\mathbb{I}_{\mathcal{A}}(\cdot)$ is the indicator function,
$\sigma^{2}=0.75$ is the variance and $k_{1}$ and $k_{2}$ are exponential
covariance kernels, each with their own unique length-scale pair. We generate
$1000$ such samples and the relative root mean square error, $\mathcal{E}$
between the DNN and FVM solutions is $3.8\%$. Fig. 20 shows a comparison of
the SPDE solutions obtained from the FVM solver and the DNN predictions for
$4$ randomly chosen stratified samples. Fig. 21 shows the histograms of
relative $L_{2}$ errors and $R^{2}$ scores of all the stratified samples.
These figures clearly show that the DNN trained on multiple length-scales data
generalizes very well to out-of-distribution input stratified data and
captures jumps at the stratification interface.

Figure 20: (Generalizability test - Case C) Each row corresponds to a randomly
chosen realization of log-input field (left column) from stratified fields
dataset and the corresponding solution response from FVM and DNN trained on
multiple length-scales GRF dataset (middle and right columns).

(a)

(b)

Figure 21: (Generalizability test - Case C) 21(a) and 21(b) corresponds to
histograms of relative $L_{2}$ errors and $R^{2}$ scores for all the input
field samples in the stratified field dataset.

###  3.5 Uncertainty Propogation (UP) and Inverse problems

The high-dimensional uncertanity propogation and inverse problems are
generally computationally expensive as a large number of forward model
evaluations of the SPDE are needed. Combining DNN approximators of SPDE
solution response with UP or inverse methods reduces this computational burden
significantly. We demonstrate this through solving few UP and inverse problems
below.

####  3.5.1 Uncertainty Propogation (UP)

Having constructed a DNN solver that estimates the solution of 1D SBVP with
good predictive accuracy in Sec. 3.2. We use this DNN solver to carry out the
following UP case-

Case 1: Fixed length-scales with $\ell_{x}=0.03$ and variance $\sigma^{2}=1$.
Next using the DNN solver trained on multiple length-scales GRF dataset (in
Sec. 3.3) we carry out three more UP cases-

Case 2: Fixed length-scales with $\ell_{x_{1}}=0.1$ and $\ell_{x_{2}}=0.3$ and
variance $\sigma^{2}=0.75$.

Case 3: Fixed length-scales with $\ell_{x_{1}}=0.06$ and $\ell_{y}=0.15$ and
variance $\sigma^{2}=0.75$.

Case 4: Uncertain length-scales with
$\ell_{x_{1}}\sim\mathrm{TN}(0.1,0.03,0.07,0.13)$ and
$\ell_{x_{2}}\sim\mathrm{TN}(0.5,0.03,0.47,0.53)$ and variance
$\sigma^{2}=0.75$, where, $\mathrm{TN}(\mu,\sigma,a,b)$ is a truncated normal
distribution with location and scale parameters, $\mu$ and $\sigma$, and
support $(a,b)$.

In each case, we draw $10^{5}$ bounded MC samples from their corresponding
input field distributions and propagate them through their respective DNN
approximators to estimate the statistics such as mean, variance and the PDF.
We compare these estimates to those obtained using the FVM solver. Fig. 22,
(23-25), (26-28) and (29-31) shows the comparison plots of output statistics
from DNN and FVM solver of all the $4$ UP cases and Tab. 4 shows the
corresponding relative $L_{2}$ error and $R^{2}$ score of the mean and
variance. Computational time taken by DNN and FVM to carry out UP is reported
in Tab. 5.

(a)

(b)

(c)

Figure 22: (UP - Case 1) 22(a), 22(b) and 22(c) corresponds to comparison
plots of mean, variance and PDF at $x=0.505$ respectively of the SBVP solution
from DNN and FVM solver for $10^{5}$ MC samples.

Figure 23: (UP - Case 2) Comparison of the mean of the SBVP solution from DNN
(left figure) and FVM solver (right figure) for $10^{5}$ MC samples.

Figure 24: (UP - Case 2) Comparison of the variance of the SBVP solution from
DNN (left figure) and FVM solver (right figure) for $10^{5}$ MC samples.

Figure 25: (UP - Case 2) Comparison of the PDF at two locations (at
$x=(0.515,0.515)$ (left figure) and $x=(0.359,0.734)$(right figure)) of the
SBVP solution from DNN and FVM solver for $10^{5}$ MC samples.

Figure 26: (UP - Case 3) Comparison of the mean of the SBVP solution from DNN
(left figure) and FVM solver (right figure) for $10^{5}$ MC samples.

Figure 27: (UP - Case 3) Comparison of the variance of the SBVP solution from
DNN (left figure) and FVM solver (right figure) for $10^{5}$ MC samples.

Figure 28: (UP - Case 3) Comparison of the PDF at two locations (at
$x=(0.515,0.515)$ (left figure) and $x=(0.359,0.734)$(right figure)) of the
SBVP solution from DNN and FVM solver for $10^{5}$ MC samples.

Figure 29: (UP - Case 4) Comparison of the mean of the SBVP solution from DNN
(left figure) and FVM solver (right figure) for $10^{5}$ MC samples.

Figure 30: (UP - Case 4) Comparison of the variance of the SBVP solution from
DNN (left figure) and FVM solver (right figure) for $10^{5}$ MC samples.

Figure 31: (UP - Case 4) Comparison of the PDF at two locations (at $x=(0.515,0.515)$ (left figure) and $x=(0.359,0.734)$(right figure)) of the SBVP solution from DNN and FVM solver for $10^{5}$ MC samples. | Mean | Variance
---|---|---
Case | Rel. $L_{2}$ error |  $R^{2}$ score | Rel. $L_{2}$ error |  $R^{2}$ score
1 | 0.00973 | 0.99929 | 0.09950 | 0.97036
2 | 0.01285 | 0.99934 | 0.08612 | 0.96895
3 | 0.01241 | 0.99938 | 0.07117 | 0.97630
4 | 0.01306 | 0.99931 | 0.09105 | 0.96519
Table 4: (UP) Relative $L_{2}$ error and $R^{2}$ score of the mean and variance of the SPDE solution for all the $4$ UP cases considered. The true statistics are estimated using $10^{5}$ MC samples of the FVM solver. Case | DNN(in hrs) | FVM(in hrs)
---|---|---
1 | 0.22 | 1.59
2 | 0.40 | 1.64
3 | 0.40 | 1.59
4 | 0.48 | 3.64
Table 5: (UP) Comparison of computational time taken by DNN and FVM solver for
UP of $10^{5}$ MC samples in all the $4$ cases.

####  3.5.2 Inverse Problem

Now we consider the problem of inferring the input field distribution based on
the observed solution response (experimental data) for a 1D SBVP in Eq. 41.

The setup for the inverse problem is as follows. We assume that the input
field, $\widetilde{A}$, is the exponential of a GRF with zero mean and an
exponential kernel whose length-scale and variance parameters are $0.03$ and
$1$ respectively. The input field is bounded between values of $0.005$ and
$33$ from below and above respectively. Similar to previous sections, let
$\Xi(\omega)$ be the flattened vector of the input field image with
$d_{\xi}=100$ pixels.

This inverse problem is solved in a Bayesian way [75, 76] as follows. $R$,
potentially noisy measurements $y_{1},y_{2},\dots,y_{R}$ of
$U(x_{1};\Xi(\omega)),\dots,U(x_{R};\Xi(\omega))$ are obtained and denoted
collectively by $D=\left\\{(x_{1},y_{1}),\dots,(x_{R},y_{R})\right\\}$.
Further we assume that these measurements are independent and identically
distributed. We model the measurement process using a Gaussian likelihood with
standard deviation $\sigma=0.032$ :

$\begin{split}{p\Big{(}D|\log(\Xi(\omega)),\theta,\sigma\Big{)}}&=\prod_{r=1}^{R}N\Big{(}y_{r}|\widehat{U}(x_{r},\Xi(\omega);\theta),\sigma^{2}\Big{)}\\\\[1.0pt]
&=\prod_{r=1}^{R}N\Big{(}y_{r}|\widehat{U}(x_{r},e^{\log(\Xi(\omega))};\theta),\sigma^{2}\Big{)},\\\\[1.0pt]
\end{split}$ (50)

where $\widehat{U}(x,\Xi(\omega);\theta)$ is our DNN trained with EF loss in
the previous 1D SBVP example (see Sec. 3.2). Our prior state-of-knowledge is
captured with a truncated normal distribution:

$p\Big{(}\log(\Xi(\omega))\Big{)}=\prod_{i=1}^{d_{\xi}}TN\Big{(}\log(\widetilde{A}(x_{i},\Xi(\omega)))|0,1,\log(0.005),\log(33)\Big{)}.$
(51)

Conditional on the observed data $D$, the posterior over the log-input field
is given by:

$p\Big{(}\log(\Xi(\omega))|D,\theta,\sigma\Big{)}=\frac{{p\Big{(}D|\log(\Xi(\omega)),\theta,\sigma\Big{)}}\cdot
p\Big{(}\log(\Xi(\omega))\Big{)}}{\mathbf{Z}},$ (52)

where $\mathbf{Z}$ is the normalizing constant. It is difficult to evaluate
$\mathbf{Z}$ analytically, but easy to sample from the likelihood and the
prior. Furthermore, the gradient of the unnormalized log posterior can be
trivially obtained. We, therefore resort to the Metropolis-adjusted Langevin
algorithm (MALA) algorithm [77, 78], a celebrated Markov-Chain Monte Carlo
(MCMC) technique which leverages the gradient of the unnormalized log
posterior to simulate a random walk whose invariant measure is the target
posterior. The step size of the random walk is set to $0.3$. The MCMC sampler
is simulated to convergence, at which point, $15,000$ random samples of
posterior are drawn by first discarding the initial 2000 samples (burn-in
phase) and thinning the rest of the sequence by 10 steps at a time to reduce
auto-correlation. We obtain an acceptance rate of $58.55\%$. Fig. 32 shows a
few samples of the input field generated from it’s posterior distribution, and
its corresponding solution responses. We see that the solutions corresponding
to the posterior sample draws of the input field capture a distribution over
the true solution, conditional on the noisy observations. The computational
time to carry out this inverse problem was around $10-15$ minutes.

Figure 32: (Inverse Problem) Green lines in the left figure and right figure
corresponds to few randomly chosen posterior samples and its corresponding
solution responses. Red lines in both left and right figures corresponds to
ground truth used to generate the experimental data. Black dots in the right
figure corresponds to experimental data/observations.

##  4 Conclusion

We developed a methodology for solving SBVPs with high-dimensional
uncertainties. The characteristic of our technique is that it does not require
a classical numerical solver (such as a finite element or finite difference
solver). Instead, we make direct use of the physics of the problem. We recast
the SBVP in a variational form and we showed that it has a unique solution in
a suitable functional space. We then minimized the corresponding energy
functional within the space of DNNs that automatically satisfy the Dirichlet
conditions. Within this functional subspace, we derived a stochastic
optimization problem with the same solution by deriving an unbiased estimator
of the energy functional. We solved this stochastic optimization problem using
a standard variant of stochastic gradient descent.

We carried out exhaustive benchmarking of our methodology by solving SBVPs
subjected to various types of high-dimensional random input fields. The
relative root mean square error, $\mathcal{E}$ between predicted and ground-
truth FVM solutions over the entirety of their particular test datasets was
between $3.86-5.3\%$ for a 2D SBVP. It was observed that the trained DNNs
generalize very well to inputs from out-of-distribution data such as
stratified fields, GRFs with lower and higher length-scales not seen in
training. In uncertainty propagation tasks, our DNNs estimated the mean with
smaller than $1.35\%$ relative $L_{2}$ error and the variance with less than
$10\%$ relative $L_{2}$ error.

Although our technique scales very well to high dimensional uncertainties in
elliptic SPDEs regardless we see an opportunity to get further improved in few
areas - First, the current computational time for training the DNNs could be
further reduced through parallelization by using multi-GPUs. Second, there is
a need for techniques that can automatically enforce the Dirichlet boundary
conditions without the need for manual construction of a trail solution to
enforce them. See for example the promising approach of [63]. Third, the
unbiased estimator we constructed to make the methodology amenable to
stochastic gradient descent was the simplest possible choice. In future work,
we will develop other estimators with smaller variance and, thus, accelerate
the convergence of stochastic optimization. Finally, and most importantly, the
behavior of the energy functional loss in the space of neural networks should
be understood and studied further analytically.

## Acknowledgements

We would like to acknowledge support from the NSF awards $\\#1737591$ and
$\\#1728165$. We would also like to acknowledge support from the Defense
Advanced Research Projects Agency (DARPA) under the Physics of Artificial
Intelligence (PAI) program (contract HR00111890034).

## References

  * Reed and Dongarra [2015] D. A. Reed, J. Dongarra,  Exascale computing and big data,  Communications of the ACM 58 (2015) 56–68.
  * Langtangen [1999] H. P. Langtangen, Computational partial differential equations: numerical methods and diffpack programming, volume 2, Springer Berlin, 1999\.
  * Smith [2013] R. C. Smith, Uncertainty quantification: theory, implementation, and applications, volume 12, Siam, 2013.
  * Sullivan [2015] T. J. Sullivan, Introduction to uncertainty quantification, volume 63, Springer, 2015\.
  * Robert and Casella [2013] C. Robert, G. Casella, Monte Carlo statistical methods, Springer Science & Business Media, 2013\.
  * Graham et al. [2018] I. G. Graham, M. J. Parkinson, R. Scheichl,  Modern monte carlo variants for uncertainty quantification in neutron transport,  in: Contemporary Computational Mathematics-A Celebration of the 80th Birthday of Ian Sloan, Springer, 2018, pp. 455–81.
  * Cliffe et al. [2011] K. A. Cliffe, M. B. Giles, R. Scheichl, A. L. Teckentrup,  Multilevel monte carlo methods and applications to elliptic pdes with random coefficients,  Computing and Visualization in Science 14 (2011) 3.
  * Kuo et al. [2012] F. Y. Kuo, C. Schwab, I. H. Sloan,  Quasi-monte carlo finite element methods for a class of elliptic partial differential equations with random coefficients,  SIAM Journal on Numerical Analysis 50 (2012) 3351–74.
  * Dick et al. [2016] J. Dick, R. N. Gantner, Q. T. L. Gia, C. Schwab,  Higher order quasi-monte carlo integration for bayesian estimation,  arXiv preprint arXiv:1602.07363 (2016).
  * Sambridge and Mosegaard [2002] M. Sambridge, K. Mosegaard,  Monte carlo methods in geophysical inverse problems,  Reviews of Geophysics 40 (2002) 3–1.
  * Sankararaman et al. [2011] S. Sankararaman, Y. Ling, S. Mahadevan,  Uncertainty quantification and model validation of fatigue crack growth prediction,  Engineering Fracture Mechanics 78 (2011) 1487–504.
  * Isukapalli et al. [1998] S. Isukapalli, A. Roy, P. Georgopoulos,  Stochastic response surface methods (srsms) for uncertainty propagation: application to environmental and biological systems,  Risk analysis 18 (1998) 351–63.
  * Angelikopoulos et al. [2012] P. Angelikopoulos, C. Papadimitriou, P. Koumoutsakos,  Bayesian uncertainty quantification and propagation in molecular dynamics simulations: a high performance computing framework,  The Journal of chemical physics 137 (2012) 144103.
  * Lockwood and Anitescu [2012] B. A. Lockwood, M. Anitescu,  Gradient-enhanced universal kriging for uncertainty propagation,  Nuclear Science and Engineering 170 (2012) 168–95.
  * Martin and Simpson [2005] J. D. Martin, T. W. Simpson,  Use of kriging models to approximate deterministic computer models,  AIAA journal 43 (2005) 853–63.
  * Bilionis and Zabaras [2012] I. Bilionis, N. Zabaras,  Multi-output local gaussian process regression: Applications to uncertainty quantification,  Journal of Computational Physics 231 (2012) 5718–46.
  * Bilionis et al. [2013] I. Bilionis, N. Zabaras, B. A. Konomi, G. Lin,  Multi-output separable gaussian process: Towards an efficient, fully bayesian paradigm for uncertainty quantification,  Journal of Computational Physics 241 (2013) 212–39.
  * Chen et al. [2015] P. Chen, N. Zabaras, I. Bilionis,  Uncertainty propagation using infinite mixture of gaussian processes and variational bayesian inference,  Journal of Computational Physics 284 (2015) 291–333.
  * Tripathy et al. [2016] R. Tripathy, I. Bilionis, M. Gonzalez,  Gaussian processes with built-in dimensionality reduction: Applications to high-dimensional uncertainty propagation,  Journal of Computational Physics 321 (2016) 191–223.
  * Najm [2009] H. N. Najm,  Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics,  Annual review of fluid mechanics 41 (2009) 35–52.
  * Eldred and Burkardt [2009] M. Eldred, J. Burkardt,  Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification,  in: 47th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, 2009, p. 976.
  * Xiu and Karniadakis [2002] D. Xiu, G. E. Karniadakis,  The wiener–askey polynomial chaos for stochastic differential equations,  SIAM journal on scientific computing 24 (2002) 619–44.
  * Ernst et al. [2012] O. G. Ernst, A. Mugler, H.-J. Starkloff, E. Ullmann,  On the convergence of generalized polynomial chaos expansions,  ESAIM: Mathematical Modelling and Numerical Analysis 46 (2012) 317–39.
  * Regis and Shoemaker [2013] R. G. Regis, C. A. Shoemaker,  Combining radial basis function surrogates and dynamic coordinate search in high-dimensional expensive black-box optimization,  Engineering Optimization 45 (2013) 529–55.
  * Volpi et al. [2015] S. Volpi, M. Diez, N. J. Gaul, H. Song, U. Iemma, K. Choi, E. F. Campana, F. Stern,  Development and validation of a dynamic metamodel based on stochastic radial basis functions and uncertainty quantification,  Structural and Multidisciplinary Optimization 51 (2015) 347–68.
  * Keogh and Mueen [2011] E. Keogh, A. Mueen,  Curse of dimensionality,  in: Encyclopedia of machine learning, Springer, 2011, pp. 257–8.
  * Constantine [2016] P. Constantine, Active subspaces: Emerging ideas for dimension reduction in parameter studies, https://speakerdeck.com/paulcon/active-subspaces-emerging-ideas-for-dimension-reduction-in-parameter-studies-1, 2016\.
  * Saltelli et al. [2000] A. Saltelli, K. Chan, E. M. Scott, et al., Sensitivity analysis, volume 1, Wiley New York, 2000\.
  * Neal [1998] R. M. Neal,  Assessing relevance determination methods using delve,  Nato Asi Series F Computer And Systems Sciences 168 (1998) 97–132.
  * Ghanem [1999] R. Ghanem,  Stochastic finite elements with multiple random non-gaussian properties,  Journal of Engineering Mechanics 125 (1999) 26–40.
  * Jolliffe [2011] I. Jolliffe,  Principal component analysis,  in: International encyclopedia of statistical science, Springer, 2011, pp. 1094–6.
  * Schölkopf et al. [1997] B. Schölkopf, A. Smola, K.-R. Müller,  Kernel principal component analysis,  in: International Conference on Artificial Neural Networks, Springer, 1997, pp. 583–8.
  * Ma and Zabaras [2011] X. Ma, N. Zabaras,  Kernel principal component analysis for stochastic input model generation,  Journal of Computational Physics 230 (2011) 7311–31.
  * Constantine et al. [2014] P. G. Constantine, E. Dow, Q. Wang,  Active subspace methods in theory and practice: applications to kriging surfaces,  SIAM Journal on Scientific Computing 36 (2014) A1500–24.
  * Constantine and Gleich [2014] P. Constantine, D. Gleich,  Computing active subspaces with monte carlo,  arXiv preprint arXiv:1408.0545 (2014).
  * Lukaczyk et al. [2014] T. W. Lukaczyk, P. Constantine, F. Palacios, J. J. Alonso,  Active subspaces for shape optimization,  in: 10th AIAA Multidisciplinary Design Optimization Conference, 2014, p. 1171.
  * Jefferson et al. [2015] J. L. Jefferson, J. M. Gilbert, P. G. Constantine, R. M. Maxwell,  Active subspaces for sensitivity analysis and dimension reduction of an integrated hydrologic model,  Computers & Geosciences 83 (2015) 127–38.
  * Constantine et al. [2015] P. G. Constantine, M. Emory, J. Larsson, G. Iaccarino,  Exploiting active subspaces to quantify uncertainty in the numerical simulation of the hyshot ii scramjet,  Journal of Computational Physics 302 (2015) 1–20.
  * Constantine et al. [2016] P. G. Constantine, C. Kent, T. Bui-Thanh,  Accelerating markov chain monte carlo with active subspaces,  SIAM Journal on Scientific Computing 38 (2016) A2779–805.
  * Tezzele et al. [2018] M. Tezzele, F. Ballarin, G. Rozza,  Combined parameter and model reduction of cardiovascular problems by means of active subspaces and pod-galerkin methods,  in: Mathematical and Numerical Modeling of the Cardiovascular System and Applications, Springer, 2018, pp. 185–207.
  * Tripathy and Bilionis [2018] R. Tripathy, I. Bilionis,  Deep uq: Learning deep neural network surrogate models for high dimensional uncertainty quantification,  arXiv preprint arXiv:1802.00850 (2018).
  * Zhu and Zabaras [2018] Y. Zhu, N. Zabaras,  Bayesian deep convolutional encoder–decoder networks for surrogate modeling and uncertainty quantification,  Journal of Computational Physics 366 (2018) 415–47.
  * Mo et al. [2018] S. Mo, Y. Zhu, J. Zabaras, Nicholas, X. Shi, J. Wu,  Deep convolutional encoder-decoder networks for uncertainty quantification of dynamic multiphase flow in heterogeneous media,  Water Resources Research (2018).
  * He et al. [2016] K. He, X. Zhang, S. Ren, J. Sun,  Deep residual learning for image recognition,  in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–8.
  * Baydin et al. [2015] A. G. Baydin, B. A. Pearlmutter, A. A. Radul,  Automatic differentiation in machine learning: a survey,  CoRR abs/1502.05767 (2015).
  * Qin et al. [2018] T. Qin, K. Wu, D. Xiu,  Data driven governing equations approximation using deep neural networks,  arXiv preprint arXiv:1811.05537 (2018).
  * Goodfellow et al. [2016] I. Goodfellow, Y. Bengio, A. Courville, Y. Bengio, Deep learning, volume 1, MIT press Cambridge, 2016.
  * Hornik et al. [1989] K. Hornik, M. Stinchcombe, H. White,  Multilayer feedforward networks are universal approximators,  Neural networks 2 (1989) 359–66.
  * Abadi et al. [2016] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al.,  Tensorflow: a system for large-scale machine learning.,  in: OSDI, volume 16, 2016, pp. 265–83.
  * Paszke et al. [2017] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, A. Lerer,  Automatic differentiation in pytorch (2017).
  * Chen et al. [2015] T. Chen, M. Li, Y. Li, M. Lin, N. Wang, M. Wang, T. Xiao, B. Xu, C. Zhang, Z. Zhang,  Mxnet: A flexible and efficient machine learning library for heterogeneous distributed systems,  arXiv preprint arXiv:1512.01274 (2015).
  * Kingma and Ba [2014] D. P. Kingma, J. Ba,  Adam: A method for stochastic optimization,  arXiv preprint arXiv:1412.6980 (2014).
  * Tieleman and Hinton [2012] T. Tieleman, G. Hinton,  Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude,  COURSERA: Neural networks for machine learning 4 (2012) 26–31.
  * Zeiler [2012] M. D. Zeiler,  Adadelta: an adaptive learning rate method,  arXiv preprint arXiv:1212.5701 (2012).
  * Raissi et al. [2017a] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations,  arXiv preprint arXiv:1711.10561 (2017a).
  * Raissi et al. [2017b] M. Raissi, P. Perdikaris, G. E. Karniadakis,  Physics informed deep learning (part ii): data-driven discovery of nonlinear partial differential equations,  arXiv preprint arXiv:1711.10566 (2017b).
  * Weinan and Yu [2018] E. Weinan, B. Yu,  The deep ritz method: A deep learning-based numerical algorithm for solving variational problems,  Communications in Mathematics and Statistics 6 (2018) 1–12.
  * Nabian and Meidani [2019] M. A. Nabian, H. Meidani,  A deep learning solution approach for high-dimensional random differential equations,  Probabilistic Engineering Mechanics 57 (2019) 14 – 25.
  * Zhu et al. [2019] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, P. Perdikaris,  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data,  arXiv preprint arXiv:1901.06314 (2019).
  * Adams and Fournier [2003] R. A. Adams, J. J. Fournier, Sobolev spaces, volume 140, Elsevier, 2003.
  * Fletcher [2013] R. Fletcher, Practical methods of optimization, John Wiley & Sons, 2013.
  * Lagaris et al. [1997] I. E. Lagaris, A. Likas, D. I. Fotiadis,  Artificial Neural Networks for Solving Ordinary and Partial Differential Equations,  ArXiv Physics e-prints (1997).
  * Berg and Nyström [2018] J. Berg, K. Nyström,  A unified deep artificial neural network approach to partial differential equations in complex geometries,  Neurocomputing 317 (2018) 28–41.
  * Szegedy et al. [2017] C. Szegedy, S. Ioffe, V. Vanhoucke, A. A. Alemi,  Inception-v4, inception-resnet and the impact of residual connections on learning.,  in: AAAI, volume 4, 2017, p. 12.
  * Wu et al. [2019] Z. Wu, C. Shen, A. Van Den Hengel,  Wider or deeper: Revisiting the resnet model for visual recognition,  Pattern Recognition (2019).
  * Veit et al. [2016] A. Veit, M. J. Wilber, S. Belongie,  Residual networks behave like ensembles of relatively shallow networks,  in: Advances in Neural Information Processing Systems, 2016, pp. 550–8.
  * Ramachandran et al. [2017] P. Ramachandran, B. Zoph, Q. V. Le,  Swish: a self-gated activation function,  arXiv preprint arXiv:1710.05941 (2017).
  * Bottou [2010] L. Bottou,  Large-scale machine learning with stochastic gradient descent,  in: Proceedings of COMPSTAT’2010, Springer, 2010, pp. 177–86.
  * Abadi et al. [2016] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. J. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Józefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. G. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. A. Tucker, V. Vanhoucke, V. Vasudevan, F. B. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng,  Tensorflow: Large-scale machine learning on heterogeneous distributed systems,  CoRR abs/1603.04467 (2016).
  * Chauvin and Rumelhart [1995] Y. Chauvin, D. E. Rumelhart, Backpropagation: theory, architectures, and applications, Psychology Press, 1995.
  * Chollet et al. [2015] F. Chollet, et al., Keras, https://keras.io/, 2015.
  * Abadi et al. [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. URL: http://tensorflow.org/, software available from tensorflow.org.
  * Guyer et al. [2009] J. E. Guyer, D. Wheeler, J. A. Warren,  Fipy: partial differential equations with python,  Computing in Science & Engineering 11 (2009).
  * Laloy et al. [2018] E. Laloy, R. Hérault, D. Jacques, N. Linde,  Training-image based geostatistical inversion using a spatial generative adversarial neural network,  Water Resources Research 54 (2018) 381–406.
  * Kaipio and Somersalo [2006] J. Kaipio, E. Somersalo, Statistical and computational inverse problems, volume 160, Springer Science & Business Media, 2006.
  * Tarantola [2005] A. Tarantola, Inverse problem theory and methods for model parameter estimation, volume 89, siam, 2005.
  * Roberts and Rosenthal [1998] G. O. Roberts, J. S. Rosenthal,  Optimal scaling of discrete approximations to langevin diffusions,  1998\.
  * Xifara et al. [2013] T. Xifara, C. Sherlock, S. Livingstone, S. Byrne, M. Girolami,  Langevin diffusions and the Metropolis-adjusted Langevin algorithm,  arXiv e-prints (2013) arXiv:1309.2983.