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

TL;DR

This paper introduces a novel deep neural network method that solves high-dimensional stochastic elliptic PDEs without relying on traditional forward solvers, effectively addressing the curse of dimensionality in uncertainty propagation.

Contribution

The authors develop a physics-informed deep residual network approach that bypasses the need for deterministic solvers, enabling efficient high-dimensional uncertainty analysis of elliptic SPDEs.

Findings

Successfully applied to various high-dimensional uncertainty scenarios

Achieved accurate uncertainty propagation without traditional forward model evaluations

Demonstrated effectiveness in inverse problems involving high-dimensional inputs

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 \textit{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 the traditional squared residual (SR) based loss function for training the ResNet, we introduce a novel 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.

Tables

Table 1. 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

Table 2. Table 2 : (1D SBVP) ResNet architecture and parameters used for each length-scale in the comparison study.

${\ell }_x$	$K$	$L$	$n$	$learningrate$	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 3. 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.

$Datasets$	$K$	$L$	$n$	Number of test samples	$ℰ$	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 4. Table 4 : (UP) Relative L 2 subscript 𝐿 2 L_{2} error and R 2 superscript 𝑅 2 R^{2} score of the mean and variance of the SPDE solution for all the 4 4 4 UP cases considered. The true statistics are estimated using 10 5 superscript 10 5 10^{5} MC samples of the FVM solver.

	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 5. Table 5 : (UP) Comparison of computational time taken by DNN and FVM solver for UP of 10 5 superscript 10 5 10^{5} MC samples in all the 4 4 4 cases.

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

Equations

$$-\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},$$

$$U(x,\omega )=G_D(x,\omega )\text{on}x\in {\Gamma }_D,$$

$$n^T(x)\left(A(x,\omega )\nabla U(x,\omega )\right)=G_N(x,\omega ),\text{on}x\in {\Gamma }_N,$$

$$\sum _{i,j=1}^d{A}_{ij}(x,\omega )v_i{v}_j\ge \alpha \parallel v{\parallel }^2,\text{for all}v\in {\mathbb{R}}^d\text{and}x\in \mathcal{X}\text{a.s.},$$

$$\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}$$

$$J[U]:=\mathbb{E}\left[{\int }_{\mathcal{X}}\left\{\frac{1}{2}\left(A\nabla U\cdot \nabla U+C{U}^2\right)-FU\right\}dx-{\int }_{{\Gamma }_N}{G}_{N}Ud{\Gamma }_N\right].$$

$$\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}$$

$$⟨U,V{⟩}_{\mathcal{H}}:=\mathbb{E}\left[{\int }_{\mathcal{X}}\left(\nabla U\cdot \nabla V+UV\right)dx\right],$$

$$\parallel U{\parallel }_{\mathcal{H}}:=⟨U,U{⟩}_{\mathcal{H}}^{\frac{1}{2}}.$$

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

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

$$T[U]=\mathbb{E}\left[{\int }_{\mathcal{X}}FUdx+{\int }_{{\Gamma }_N}{G}_{N}Ud{\Gamma }_N\right],$$

$$\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}$$

$$\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}$$

$$B[V,U]=T[V],$$

$$\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}$$

$$B[v\Phi ,U^{\ast }]$$

$$\Rightarrow \mathbb{E}\left[{\int }_{\mathcal{X}}\left(A\nabla U^{\ast }\cdot \nabla (v\Phi )+C{U}^{\ast }(v\Phi )\right)dx\right]$$

$$\Rightarrow \mathbb{E}\left[\Phi {\int }_{\mathcal{X}}\left(A\nabla U^{\ast }\cdot \nabla v+C{U}^{\ast }v\right)dx\right]$$

$${\int }_{\mathcal{X}}\left(A\nabla U^{\ast }\cdot \nabla v+C{U}^{\ast }v\right)dx={\int }_{\mathcal{X}}Fvdx+{\int }_{{\Gamma }_N}{G}_{N}vd{\Gamma }_N\text{a.s.},$$

$$\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}$$

$$\begin{array}{r}-\nabla \cdot \left(\tilde{A}(x,\omega )\nabla U(x,\omega )\right)+C(x)U(x,\omega )=F(x),\text{in}x\in \mathcal{X},\\ U(x,\omega )=G_D(x)\text{on}x\in {\Gamma }_D,\\ n^T(x)\left(\tilde{A}(x,\omega )\nabla U(x,\omega )\right)=G_N(x),\text{on}x\in {\Gamma }_N.\end{array}$$

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

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

$$z^{(1,0)}=W_0^{(1,0)}x+W_1^{(1,0)}\Xi (\omega )+b^{(1,0)},$$

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

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

$$z^{(i,0)}:=z^{(i,L)}+z^{(i-1,0)},$$

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

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

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

Full text

††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.

[FIGURE:]

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):

[TABLE]

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

[TABLE]

a.s. and

[TABLE]

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

[TABLE]

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

[TABLE]

and this solution minimizes the functional:

[TABLE]

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

[TABLE]

with inner product

[TABLE]

and the corresponding norm

[TABLE]

Notice that

[TABLE]

where

[TABLE]

is a bounded bilinear form and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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,

[TABLE]

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],

[TABLE]

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].

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:

[TABLE]

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:

[TABLE]

followed by

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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 :

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

[TABLE]

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 -

[TABLE]

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

[TABLE]

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:

[TABLE]

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))

[TABLE]

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.

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.

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.

3.3 Stochastic boundary value problem in 2D

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

[TABLE]

with boundary conditions:

[TABLE]

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.

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

the mean and covariance functions are chosen as follows:

[TABLE]

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.

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:

[TABLE]

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.

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.

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.

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.

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.

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:

[TABLE]

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.

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.

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$ :

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

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).