Deep convolutional neural networks for uncertainty propagation in random fields

TL;DR

This paper introduces a hierarchical deep convolutional neural network surrogate model for high-dimensional uncertainty propagation, effectively capturing complex input-output relationships in continuum mechanics problems.

Contribution

The study presents a novel hierarchical CNN architecture tailored for high-dimensional uncertainty quantification, enhancing efficiency and accuracy over traditional methods.

Findings

The surrogate accurately models diverse I/O mappings.

It effectively characterizes statistical properties of output fields.

Numerical results outperform Monte Carlo estimates.

Abstract

The development of a reliable and robust surrogate model is often constrained by the dimensionality of the problem. For a system with high-dimensional inputs/outputs (I/O), conventional approaches usually use a low-dimensional manifold to describe the high-dimensional system, where the I/O data is first reduced to more manageable dimensions and then the condensed representation is used for surrogate modeling. In this study, we present a new solution scheme for this type of problems based on a deep learning approach. The proposed surrogate is based on a particular network architecture, i.e. the convolutional neural networks. The surrogate architecture is designed in a hierarchical style containing three different levels of model structures, advancing the efficiency and effectiveness of the model in the aspect of training and deploying. To assess the model performance, we carry out uncertainty quantification in a continuum mechanics benchmark problem. Numerical results suggest the proposed model is capable of directly inferring a wide variety of I/O mapping relationships. Uncertainty analysis results obtained via the proposed surrogate have successfully characterized the statistical properties of the output fields compared to the Monte Carlo estimates.

Figures

Tables

Table 1. Table 1: Summary of the case study analysis

Case study	Inputs	Outputs	Surrogate modeling
Case $1$	$𝒙=E\left(s\right)$	$𝒚=w\left(s\right)$	$\hat{f}(\cdot ):{ℝ}^{4096}\to {ℝ}^{4096}$
Case $2$	$𝒙=E\left(s\right)$	$𝒚={[{\sigma }_v\left(s\right),{\tau }_{\max }\left(s\right),{\tau }_{xy}\left(s\right)]}^{†}$	$\hat{f}(\cdot ):{ℝ}^{4096}\to {ℝ}^{12288}$
Case $3$	$𝒙=[E\left(s\right),f\left(s\right)]$	$𝒚=[w\left(s\right),{\sigma }_v\left(s\right)]$	$\hat{f}(\cdot ):{ℝ}^{8192}\to {ℝ}^{8192}$

Equations

$${\mathcal{L}}_{\mathbf{s}}\left(\mathbf{\Theta }(\mathbf{s}),\mathbf{u}(\mathbf{s})\right)=\mathbf{f}(\mathbf{s})\text{in}\mathcal{S}$$

$$\mathbf{u}=g^{(D)}$$

$$\nabla \mathbf{u}\cdot \hat{n}=g^{(N)}$$

$$\mathbf{u}:\mathcal{S}\times \Omega ⟶{\mathbb{R}}^d$$

$${\mathbf{\mu }}_f={\int }_{\Omega }f(\mathbf{x})p(\mathbf{x})\mathrm{d}\mathbf{x}$$

$${\mathbf{\sigma }}_f^2={\int }_{\Omega }{\left(f(\mathbf{x})-{\mathbf{\mu }}_f\right)}^{2}p(\mathbf{x})\mathrm{d}\mathbf{x}$$

$$p_f(\mathbf{y})={\int }_{\Omega }\delta (\mathbf{y}-f(\mathbf{x}))p(\mathbf{x})\mathrm{d}\mathbf{x}$$

$${\mathbf{\mu }}_f$$

$${\mathbf{\sigma }}_f^2$$

$$p_f(\mathbf{y})$$

$${\mathbf{\gamma }}_{m,n}=(\mathbf{x}\star \mathbf{W}{)}_{m,n}^s=\sum _{p=1}^h\sum _{q=1}^w\sum _{r=1}^C{\mathbf{W}}_{p,q,r}^s\cdot {\mathbf{x}}_{p+m,q+n,r}$$

$$W^{\prime }=\frac{W-w+2p}{s}+1\text{and}{H}^{\prime }=\frac{H-h+2p}{s}+1$$

$$h(\boldsymbol{\gamma}_{m,n})=\max(\boldsymbol{\gamma}_{m,n},\epsilon)\quad\text{and}\quad h^{\prime}(\boldsymbol{\gamma}_{m,n})=\left\{\begin{array}[]{ll}{1}&\text{if}\quad{\boldsymbol{\gamma}_{m,n}>\epsilon}\\ {0}&\text{if}\quad{\boldsymbol{\gamma}_{m,n}\leq\epsilon}\end{array}\right.$$

$$\mathbf{\vartheta }=\frac{\alpha \left(h(\mathbf{\gamma })-{\mu }_{\mathbf{x}}\right)}{\sqrt{{\sigma }_{\mathbf{x}}^2+ϵ}}+\beta$$

$$\hat{f}(\mathbf{x})={{\rm Y}}^K\circ {{\rm Y}}^{K-1}\circ \dots {{\rm Y}}^1(\mathbf{x})$$

$${\mathbf{x}}_l={\mathbf{y}}_{l-1}={{\rm Y}}_{l-1}\left({\mathbf{x}}_{l-1}\right)$$

$${\mathbf{x}}_l={\mathbf{y}}_{l-1}={{\rm Y}}_{l-1}\left(\left[{\mathbf{x}}_0,{\mathbf{x}}_1,\dots {\mathbf{x}}_{l-1}\right]\right)$$

$$\mathcal{L}\left(\mathbf{\theta };\mathcal{D}\right)=\frac{1}{N}\sum _{i=1}^N{\Vert {\mathbf{y}}_i-\hat{f}\left({\mathbf{x}}_i,\mathbf{\theta }\right)\Vert }^2$$

$$\mathcal{L}\left(\mathbf{\theta };\mathcal{D}\right)=\frac{1}{N}\sum _{i=1}^N{\Vert {\mathbf{y}}_i-\hat{f}\left({\mathbf{x}}_i,\mathbf{\theta }\right)\Vert }^2+\lambda \sum _{j=1}^K{\Vert {\mathbf{W}}^{(j)}\Vert }_2^2$$

$${\mathbf{\theta }}_{k+1}\leftarrow {\mathbf{\theta }}_k-{\eta }_k{\nabla }_{\mathbf{\theta }}\mathcal{L}\left(\mathbf{\theta };{\mathcal{D}}_M\right)$$

$${\mathbf{\theta }}_{k+1}\leftarrow {\mathbf{\theta }}_k+{\eta }_k\frac{M_k}{1-{\beta }_1^k}/\left(\sqrt{\frac{V_k}{1-{\beta }_2^k}}+ϵ\right)$$

$$M_k$$

$$V_k$$

$$RMSE=\sqrt{\frac{1}{M}\sum _{i=1}^M||{\hat{\mathbf{y}}}_i-{\mathbf{y}}_i|{|}_2^2}$$

$$R^2=1-\frac{Var\left(\mathbf{y}|\mathbf{x}\right)}{Var\left(\mathbf{y}\right)}=1-\frac{{\sum }_{i=1}^M||{\hat{\mathbf{y}}}_i-{\mathbf{y}}_i|{|}_2^2}{{\sum }_{i=1}^M||{\bar{\mathbf{y}}}_i-{\mathbf{y}}_i|{|}_2^2}$$

$$-\mathrm{div}\mathbf{C}\epsilon (\mathbf{\theta })-\gamma =0$$

$$-\mathrm{div}\gamma =f$$

$$-\lambda t^{-2}(\nabla w-\mathbf{\theta })$$

$$\mathbf{\theta }=0,w=0$$

$$\log \mathbf{x}(s)\sim \mathrm{GP}\left(m(s),k\left(s,s^{\prime }\right)\right)$$

$$k\left(s,s^{\prime }\right)={\sigma }^2\exp \left(-\frac{{\left(s-s^{\prime }\right)}^2}{2{\ell }^2}\right)$$

$$64\times 64\times 2$$

$$\displaystyle 64\times 64\times 2=\begin{array}[]{ll}64\times 64\times 1\\ 64\times 64\times 1\end{array}$$

Full text

\draftfalse\journalname

Wiley

NatHaz Modeling Laboratory, University of Notre Dame, Notre Dame, IN, USA.

\correspondingauthor

Xihaier Luoxluo1@nd.edu

\justify

{keypoints}

A surrogate model for quantifying the effect of spatially random attributes in coupled elliptic systems is developed.

A hierarchical modeling design using deep convolutional neural networks as the underlying components is proposed to improve the model efficiency in terms of training and deploying.

The efficacy of the proposed learning approach is demonstrated using a benchmark structural problem with a variety of high-dimensional mapping relationships: $\mathbb{R}^{4096}\rightarrow\mathbb{R}^{4096},\mathbb{R}^{4096}\rightarrow\mathbb{R}^{12288},\text{and }\mathbb{R}^{8192}\rightarrow\mathbb{R}^{8192}$.

1 Introduction

The primitive parameters of real-world mechanical systems such as material properties and external excitations are often not known exactly and are subject to uncertainty. This uncertainty may arise from inevitable manufacturing imperfections, incomplete knowledge on system parametrization, and others. In most cases, uncertain parameters are modeled as random variables utilizing a sequence of probability distributions and the system of interest is characterized by stochastic partial differential equations (SPDEs) Smith (\APACyear2013). Because the input parameters are random, the output of the system also becomes a random variable. To obtain comprehensive probabilistic descriptions for the random output quantities, one can rely on the ensemble of independent samples drawn from the distribution of the inputs using the Monte Carlo method, evaluating the system at a finite number of realizations and estimating the statistics of the target outputs (See Fig. 1 .(A)). Such uncertainty quantification (UQ) task involving many system queries can be greatly expedited through the use of surrogate models especially for applications that depend on a time-demading simulator Viana \BOthers. (\APACyear2014).

Conceptually simple, a formidable challenge remains open concerning the development of a reliable and robust surrogate in cases where the input/output is no longer a random variable, instead, is a random field or even a set of random fields. For instance, one can treat material properties as random variables, which are random in essence but are uniformly distributed over the domain. Such treatment ignores the possible spatial variation and the realization of perfectly uniform material is unlikely to occur in practice, thus fails to provide a realistic reflection of the real-world situations Babuška \BOthers. (\APACyear2007); Gurley \BBA Kareem (\APACyear1999); G. Chen \BOthers. (\APACyear2016). On the other hand, the random field theory offers a more reliable avenue as it systematically takes the spatial correlation into accounts (See Fig. 1 .(B)). The numerical discretization of a random field with small correlation can easily result in a large number of random variables. Unfortunately, most existing surrogate models dealing with uncertainty propagation tasks have difficulty scaling to high-dimensional problems and the conjunction use of the random field and the time-demanding computer model can be prohibitive Schwab \BBA Gittelson (\APACyear2011); Babuška \BOthers. (\APACyear2007); Dai \BBA Cao (\APACyear2017).

Given that the direct construction of surrogate models yielding an easy-to-evaluate mechanism between high-dimensional inputs and outputs can be extremely demanding, practitioners usually attempt to project the given random field data from its original high dimension to a lower dimension. The standard approach to performing dimensionality reduction on the continuous functional space in the UQ literature is the Karhunen-Loève expansion (KLE) Xiu \BBA Karniadakis (\APACyear2002); Ghanem \BBA Spanos (\APACyear2003); Le Maître \BBA Knio (\APACyear2010). Known under different names such as principal component analysis, proper orthogonal decomposition, etc, the KLE first projects a random field, which is infinite dimensional, to finite dimensions using eigenfunctions of the random field’s covariance function as basis functions. Then, the KLE represents the random field as a linear combination of these orthogonal eigenfunctions and their corresponding uncorrelated coefficients, which are random variables. In spite of numerous applications in which the KLE based dimensionality reduction scheme circumvents the curse of dimensionality, the KLE by definition can only preserve the first and second order statistics, i.e., the mean and covariance of the random field. However, uncertainties arise from material heterogeneities or randomly distributed force fields usually exhibit nonlinear structures Kareem (\APACyear2008). Moreover, the KLE is an unsupervised learning algorithm from a machine learning perspective Murphy (\APACyear2012), indicating it does not explicitly integrate the output information into the dimensionality reduction procedure and additional classification/regression model is needed to complete the surrogate modeling task (See the upper green route in Fig. 1. (C)).

In recent years, deep neural networks (DNNs) have gained increased popularity in representing the high-dimensional solutions of a wide class of partial differential equations T\BPBIQ. Chen \BOthers. (\APACyear2018); Raissi \BOthers. (\APACyear2019); Sirignano \BBA Spiliopoulos (\APACyear2018); Mo \BOthers. (\APACyear2019); Zhu \BBA Zabaras (\APACyear2018). Compared to the conventional surrogate modeling approaches, DNNs tackle the curse of dimensionality by learning the latent representation through the use of a series of layers, which are connected by affine operations and nonlinear activations. Usually formed as a supervised learning problem, DNNs are capable of providing an end-to-end field-to-field surrogate model, thus avoiding the need to build a dimension reduction model (See the lower red route in Fig. 1. (C)). In this work, we attempt to deploy one of the most widely used network architectures, the convolutional neural networks (CNNs), to the UQ analysis of civil structures with spatially varying system properties. Within the realm of neural networks modeling, CNNs are the workhorse of image classification Gao \BBA Mosalam (\APACyear2018); Krizhevsky \BOthers. (\APACyear2012), object detection Cha \BOthers. (\APACyear2017); Gopalakrishnan \BOthers. (\APACyear2017), computer vision Yeum \BBA Dyke (\APACyear2015); Yang \BOthers. (\APACyear2018), and signal processing Rafiei \BBA Adeli (\APACyear2018). This type of neural networks gains computational efficiency in dealing with data of a grid-like topology by extracting multi-scale features within the modeling process. This unique learning property triggers the use of CNNs for surrogate modeling of the complex nonlinear mapping between discretized random fields.

In this paper, a machine learning based structural analysis framework is proposed to statistically characterize the system responses under various spatially varying properties. We formulate our surrogate model in a hierichiacal form where CNNs are adopted as fundamental model components. Because CNNs are prone to overparameterization, additional network connections between nonadjacent convolutional layers are created in such a way that extracted features can be reused, thus reducing the model parameters and strengthening the machine learning process He \BOthers. (\APACyear2016); Huang \BOthers. (\APACyear2017). Moreover, the operations of deconvolution included in the CNNs can generate checkerboard artifacts, especially in the context of regression problems of high-resolution images. To address this issue, we make the joint use of the upsampling methods such as bicubic resizing and convolution operation in the decoding process instead of conventional deconvolution Odena \BOthers. (\APACyear2016); Dong \BOthers. (\APACyear2015). To make the training process more efficient, we incorporate the annealing learning strategy in the setting of stochastic gradient descent and adopt batch normalization method to avoid data distortion Robbins \BBA Monro (\APACyear1951); Ioffe \BBA Szegedy (\APACyear2015). The overall machine learning methodology is applied to a plate-like structure governed by a set of fourth-order elliptic partial differential equations. We demonstrate the flexibility and robustness of our proposed surrogate model by considering one-to-one ($\mathbb{R}^{4096}\rightarrow\mathbb{R}^{4096}$), one-to-many ($\mathbb{R}^{4096}\rightarrow\mathbb{R}^{12288}$), and many-to-many ($\mathbb{R}^{8192}\rightarrow\mathbb{R}^{8192}$) mappings, where each one represents a discretized high-dimensional random field.

In the rest of this paper, Section 2 gives the problem statement for using deep neural networks to propagate uncertainties represented by random fields. Section 3 provides a comprehensive guide in terms of building, training, and validating the proposed surrogate. In Section 4, a systematic case study on first identifying intrinsic relationships between discretized random fields and then applying the learned surrogate to the uncertainty analysis is presented. Lastly, Section 5 summarizes the principal conclusions and discusses the future work.

2 Problem statement

This paper is focused on the development of surrogate models for effectively propagating and quantifying the influence of spatially varying uncertainties on structural behaviors in applied mechanics problems. To fulfill this objective, we consider general boundary value problems of the following elliptic PDEs form:

[TABLE]

where $\boldsymbol{s}$ is the cartesian coordinates, $\mathcal{L}_{\boldsymbol{s}}$ is a differential operator, $\boldsymbol{\Theta}$ represents an input property field encompassed in the constitutive equation, $\boldsymbol{u}$ denotes an unknown solution to the equation, $\boldsymbol{f}$ represents the source term that drives the system, and $\mathcal{S}$ is a Lipschitz smooth domain. To determine a unique solution, boundary conditions, i.e. the Dirichlet and Neumann boundary conditions are imposed on $\Gamma_{\mathrm{D}}\subset\partial\mathcal{S}$ and $\Gamma_{\mathrm{N}}\subset\partial\mathcal{S}$, respectively:

[TABLE]

with $\hat{n}$ denoting the unit normal on the boundary. The operator $\nabla$ is considered only with respect to (w.r.t.) the spatial variables $\boldsymbol{s}$. In a typical probabilistic mechanics problem, one often deals with systems subjected to stochastic excitations $\boldsymbol{f}(\boldsymbol{s})$, or themselves involve random parameters $\boldsymbol{\Theta}(\boldsymbol{s})$, or both. A standard uncertainty analysis framework for quantifying the influence of system uncertainties comprise three parts, that is, uncertainty representation, uncertainty propagation, and uncertainty quantification Smith (\APACyear2013).

Remark* (1).*

Without loss of generality, we formulate the problem (Eq. 1) using a general elliptic operator $\mathcal{L}_{\boldsymbol{s}}$. In the context of uncertainty quantification of structural systems, the operator $\mathcal{L}_{\boldsymbol{s}}$ could be the momentum balance equation, which reads the mass density, body force, initial velocity, etc as the inputs and computes the corresponding displacement field. And our goal is to characterize and manage the spatial variability induced uncertainty effects on the system performance.

2.1 Step 1: random fields for uncertainty representation

Let $(\Omega,\mathscr{S},\mathbb{P})$ be a complete probability space, where $\Omega$ is the set of the elementary events $\omega$, $\mathscr{S}$ is an $\sigma$-algebra on the set $\Omega$, and $\mathbb{P}$ is a probability measure. We are interested in quantifying the uncertainty associated with the solution $\boldsymbol{u}$ that is due to the spatial variability of the inputs $\boldsymbol{\Theta},\boldsymbol{f}$. To achieve the goal, we use random fields, which are second order stationary and statistically homogeneous, to represent physically relevant variability Shinozuka \BBA Jan (\APACyear1972). In general, a random field $\boldsymbol{r}(\boldsymbol{s},\omega)$ is a manifold in the Hilbert space consisting of a family of random variables. Hence, the solution mapping of Eq. 1 in a geometrically bounded regime with discretized dimension $d$ is expressed as:

[TABLE]

The central work of the first uncertainty representation step is to generate random field samples. Usually, the random field generation covers two parts, where the first part focuses on generating a set of uncorrelated random variables and the second part transfers these random variables into a correlated random field with prescribed statistical properties. There are several methods available for this purpose and matrix decomposition method is considered in this paper Davis (\APACyear1987); Constantine \BBA Wang (\APACyear2010); Kareem (\APACyear2008). Specifically, the realization of a random field with known properties is obtained by multiplying the lower triangular matrix computed from the Cholesky decomposition of the prescribed covariance matrix of $\boldsymbol{r}(\boldsymbol{s},\omega)$ by a random vector comprising a set of standard normal random variables.

2.2 Step 2: deep neural networks for uncertainty propagation

Given the fact that the input uncertainty is represented by means of random fields, the second step of the uncertainty analysis is to propagate the spatial variability to the outputs. To achieve this goal, deep neural networks are considered in this paper for two main reasons. First, the performance function $f\left(\cdot\right)$ of structural systems with spatially varying properties is highly complex and implicit in general. To approximate the function $f\left(\cdot\right)$, neural networks based surrogate is a strong candidate since it has been shown that a single hidden layer based neural networks are universal function approximators under mild conditions Hornik \BOthers. (\APACyear1989). Secondly, going deep allows the network structure smoothly reduces the high dimensionality $d$ arises from the discretization of the random field Bengio (\APACyear2012).

Mathematically, the approach of using deep neural networks to propagate random fields represented uncertainties can be summarized as follows: (1) collect the training dataset $\mathcal{D}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{n}$ by querying $f\left(\cdot\right)$ with $\boldsymbol{x}$ and $\boldsymbol{y}$ denoting the input and output respectively. For instance, $\boldsymbol{x}$ can be the system property field $\boldsymbol{\Theta}(\boldsymbol{s})$ or the external loading field $\boldsymbol{f}(\boldsymbol{s})$ or both $\boldsymbol{x}=[\boldsymbol{\Theta},\boldsymbol{f}]$ , and $\boldsymbol{y}$ usually refers to the solution field $\boldsymbol{u}$. (2) train the surrogate model using $\mathcal{D}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{n}$. The central goal is to establish a functional relationship $\hat{f}\left(\cdot\right)$ in such a way that discrepancies between $f\left(\cdot\right)$ and $\hat{f}\left(\cdot\right)$ are optimally minimized. Once $\hat{f}\left(\cdot\right)$ is available, one can immediately use it for the uncertainty quantification.

2.3 Step 3: Monte Carlo method for uncertainty quantification

With the aid of an efficient surrogate model, uncertainties concerning the inputs can be effectively propagated into the system outputs. Subsequently, the last uncertainty quantification step centers on estimating the statistical properties of the output. The most widely used objective of a standard UQ problem is the characterization of the first two statistical moments, i.e., the mean,

[TABLE]

and the variance,

[TABLE]

In this paper, we are also interested in the computation of the probability density function (PDF):

[TABLE]

Since the function $f\left(\cdot\right)$ may not always be known in closed form for the problem under consideration, the sampling methods such as the Monte Carlo (MC) method are used to numerical approximations of the multidimensional integrals Liu (\APACyear2008); Robert \BBA Casella (\APACyear2013). In particular, the high fidelity model $f\left(\cdot\right)$ is replaced by the deep neural networks based surrogate model $\hat{f}\left(\cdot\right)$. And the integral for the statistical moments as well as the PDF can be consequently expressed in the form of the averaging operator using sampled values:

[TABLE]

The noteworthy fact is the combined use of the surrogate model and MC method in the context of the uncertainty quantification can usually reduce the computational burden by orders of magnitude. More importantly, the construction of a reliable and robust surrogate model is of critical significance throughout the uncertainty analysis. In the following section, we propose a deep neural networks based surrogate to approximate the mapping relationship between discretized high-dimensional random fields.

3 Methodology

The fundamental structure of the proposed surrogate model is based on a specific neural networks architecture, that is, the convolutional neural networks (CNNs) LeCun \BOthers. (\APACyear1998). We build our surrogate in a hierarchical manner so that the model is more efficient in terms of training and deploying. Because we mainly focus on the regression problems of random fields, we name the model field regressor. The rest of this section covers a guideline on the construction, training, and validation of a field regressor.

3.1 Model construction

Given the input data set $\boldsymbol{x}$ and the desired output data set $\boldsymbol{y}$, we aim to build a machine-learned model $\hat{f}\left(\cdot\right)$ such that $\hat{f}\left(\cdot\right)$ can approximate the input-output relationship $\boldsymbol{y}=f\left(\boldsymbol{x}\right)$. In the present case, $\boldsymbol{x}$ and $\boldsymbol{y}$ are discretized random fields. Let $\boldsymbol{x}\in\mathbb{R}^{H\times W\times C}$ and $\boldsymbol{y}\in\mathbb{R}^{H^{\prime}\times W^{\prime}\times C^{\prime}}$, where $H(H^{\prime})$ is the y-axis resolution, $W(W^{\prime})$ is the x-axis resolution, and $C(C^{\prime})$ denotes the number of random fields. The proposed surrogate $\hat{f}\left(\cdot\right)$ is essentially a hierarchical model that contains three different levels. It approximates the nonlinear regression function $f\left(\cdot\right)$ through a hierarchy of convolutional layers. Fig. 2 summarizes the specifications for the architecture and assembling of the surrogate $\hat{f}\left(\cdot\right)$. We will investigate the detailed designs in a bottom-up manner, starting from level $1$ model.

3.1.1 Level 1 model: convolutional layer

A convolutional layer by its most prevailing definition is a composite function that has three mathematical operations Dumoulin \BBA Visin (\APACyear2016). First, a convolutional layer performs element-wise multiplications between a subarray of the input field $\boldsymbol{x}$ and the kernel matrix $\mathbf{W}\in\mathbb{R}^{h\times w\times C}$ to obtain the so-called feature value at location $(m,n)$:

[TABLE]

where $s$ denotes the kernel number, i.e. $s=1,2,\dots,C^{\prime}$. The computed output matrices are often called the feature maps, which contain the local correlation information about the input field. There are two other important parameters for the convolutional layer: the stride number $s$ and the padding number $p$. Specifically, $s$ represents the distance between two consecutive kernels and $p$ introduces zeros values to the borders of the input fields. The selection of $h$, $w$, $s$, and $p$ together determines the size of the output fields Goodfellow \BOthers. (\APACyear2016):

[TABLE]

Fig. 3 gives an example illustrating the convolution arithmetic. In this example, the input field $\boldsymbol{x}\in\mathbb{R}^{5\times 5\times 1}$, the kernel matrix $\mathbf{W}\in\mathbb{R}^{3\times 3\times 1}$, the stride number $s=2$ and the padding number $p=1$. Hence, the output field takes the form of $\boldsymbol{y}\in\mathbb{R}^{3\times 3\times 1}$.

Secondly, the output of the convolution flows over an element-wise nonlinear function that is often referred to as the activation function Goodfellow \BOthers. (\APACyear2016). Numerous types of activation functions exist in the literature and popular choice include the Sigmoid, the hyperbolic tangent or the rectified linear unit (ReLU). In conjunction with empirical knowledge, ReLU is implemented in the proposed scheme for two reasons. It has better gradient propagation property and analytical expressions of the ReLU function and its gradient are easy-to-implement as they only require basic addition and multiplication. An element-wise expression is given as:

[TABLE]

with $\epsilon$ denoting a small number, typically [math]. By manipulating the network architecture (e.g. number of layers, size of each layer, etc.), one can learn functions of arbitrary complexity.

Thirdly, we apply the batch normalization to the output of Eq. 10:

[TABLE]

where $\alpha$ and $\beta$ are parameters to be learned, $\mu_{\boldsymbol{x}}$ and $\sigma_{\boldsymbol{x}}^{2}$ are sample mean and variance respectively. By means of Eq. 11, the internal covariate shift of the network architecture is minimized, allowing the data to propagate through a deeply structured neural network without distortion Ioffe \BBA Szegedy (\APACyear2015).

3.1.2 Level 2 model: dense block

A quick summary, the level 1 model, i.e. the aforestated convolutional layer is structured as $\Upsilon(\cdot):BN\circ ReLU\circ Conv(\boldsymbol{x})$. Neural networks approximate the input-output relationship $f(\cdot)$ as:

[TABLE]

where $K$ is the number of layers defined in the network architecture and $\circ$ is taken to be the composition operator with $\Upsilon^{j},j\in\{1,2,\cdots,K\}$. Usually, CNNs take the output of the $(l-1)^{th}$ convolutional layer as the input to the $(l)^{th}$ layer:

[TABLE]

Though simple and neat, recent research reveals such block design often encounters information degradation problems, which impede the training of CNNs Glorot \BBA Bengio (\APACyear2010). To remedy this issue, distant connections among nonadjacent layers are established, where every single layer is fully connected with all the subsequent layers He \BOthers. (\APACyear2016); Huang \BOthers. (\APACyear2017). Consequently, the improved expression of function writes as:

[TABLE]

In Fig. 2. (D), these augmented connections are depicted by the dashed lines. And a sequence of convolutional layers with nonadjacent layers is called dense block. It is assumed that each $\Upsilon_{i}(i=1,2,\ldots l-1)$ extracts $k$ feature maps. As a result, the input field for the $(l)^{th}$ layer has a size of $\mathbb{R}^{W_{l}\times H_{l}\times C_{l}}$, where $C_{l}=k_{0}+k(l-1)$. Noticeably, $k_{0}$ denotes the number of random fields of the first layer and $k$ is the growth rate.

3.1.3 Level 3 model: sampling processor

To improve the efficiency of $\hat{f}(\cdot)$ in terms of training, the high-dimensional input data $\boldsymbol{x}$ will first go through a downsampling processing, where feature maps produced by different convolutional layers are continually downsized to a coarser spatial resolution $H_{l}<H_{l-1}\,\,\&\,\,W_{l}<W_{l-1}$, resulting in a significant reduction of model parameters. This can be achieved either by controlling the stride number $s$ or using max pooling operation. In this paper, we choose to reduce the size of $\boldsymbol{x}$ by adopting non-unit stride so that sub-regional data relations can be better preserved Dumoulin \BBA Visin (\APACyear2016). On the other hand, the dimension-reduced data will be projected back to the high dimension output space. Projection methods mainly fall into two categories, i.e. the deconvolution method and the interpolation method. We choose to resize the data by means of bicubic interpolation since the deconvolution method tends to produce checkerboard artifacts Odena \BOthers. (\APACyear2016).

3.2 Model training

3.2.1 Regularized loss function

Training a field regressor centers on optimizing the network parameters $\boldsymbol{\theta}$ in such a way that the discrepancy between model predictions $\hat{\boldsymbol{y}}_{i}$ and true observations $\boldsymbol{y}_{i}$ is minimized. For regression tasks, the discrepancy is typically measured by the mean squared error (MSE) loss. Given a training dataset $\mathcal{D}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{N}$ consisting of $n$ independent and identically distributed (i.i.d.) samples, the loss function is defined as:

[TABLE]

In parctice, neural nets are notorious for overfitting. One usually resort to regularization techniques, penalizing $f^{{}^{\prime}}(\cdot)$ by either promoting sparsity or driving weights to [math]. In this paper, $L_{2}$ norm is chosen to penalize the model complexity Goodfellow \BOthers. (\APACyear2016). The regularized MSE loss writes as:

[TABLE]

where $K$ is the total number of layers, and $\lambda$ is a small constant that forces $f^{{}^{\prime}}(\cdot)$ learning small weights (often termed as the weight decay), and $\mathbf{W}^{j}$ denotes all the parameters in the $j^{th}$ layer. For the field regressor, $\mathbf{W}^{j}$ includes the kernel weights utilized in the convolution operation, the scale and shift parameters in the batch normalization operation, and bicubic interpolation parameters used in the image resampling.

3.2.2 Mini-batch stochastic gradient descent

In order to compute the parameters $\boldsymbol{\theta}$, the gradient-based optimization is adopted as $f^{{}^{\prime}}(\cdot)$ is a highly composite function involving a series of layers of nonlinear transformation Rumelhart \BOthers. (\APACyear1988). In particular, the Mini-batch stochastic gradient descent (mini-batch SGD) is implemented to solve this problem Goodfellow \BOthers. (\APACyear2016). The central idea is to approximate the negative gradient of the loss function by means of the expectation of a batch of samples:

[TABLE]

where $\eta_{k}$ is the learning rate in the $k^{th}$ iteration and $\mathcal{L}\left(\boldsymbol{\theta};\mathcal{D}_{M}\right)$ is the loss function evaluated using a mini-batch subset of samples $\mathcal{D}_{M}\subset\mathcal{D}$. Specifically, a mini-batch is initialized between $5$ and $500$ samples, chosen through a random manner and the Adaptive Moments (ADAM) optimization algorithm, a variant of the SGD methods, is selected to update the learnable parameters Kingma \BBA Ba (\APACyear2014):

[TABLE]

In Eq. 18, $\epsilon$ is a small number introduced to prevent [math] denominator, and $M_{k}$ and $V_{k}$ are estimates of the mean and variance of the gradients, respectively. In ADAM, the update scheme of them are given as follows:

[TABLE]

Note that $M_{0}$ and $V_{0}$ are set to [math] during the initialization. $\beta_{1}$ and $\beta_{2}$ are averaging parameters. In this work, values are configured as $\beta_{1}=0.9$ and $\beta_{2}=0.999$ following the suggestion in Kingma \BBA Ba (\APACyear2014).

3.3 Model validation

To assess the performance of a fully trained model, $f^{{}^{\prime}}\left(\cdot\right)$ is tested on a validation dataset $\mathcal{D}=\{\boldsymbol{x_{i}},\boldsymbol{y_{i}}\}_{i=1}^{M}$. More specifically, two performance indicators are selected to check whether $f^{{}^{\prime}}\left(\cdot\right)$ provides a good approximation. First, root mean squared error (RMSE) is used to monitor the convergence history of model training:

[TABLE]

Secondly, the coefficient of determination ($R^{2}$) is selected as the performance indicator. In statistics, $R^{2}$ estimates the proportion of variance in $\boldsymbol{y}$ that can be explained by $f^{{}^{\prime}}\left(\cdot\right)$:

[TABLE]

where $\bar{\boldsymbol{y}}_{i}$ is the mean value. From a mathematical perspective, $R^{2}$ takes on values between 0 and 1. When $R^{2}$ is close to one, it indicates $f^{{}^{\prime}}\left(\cdot\right)$ fits the data well. On the other hand, $R^{2}=0$ means the model does not explain any variability in $\boldsymbol{y}$. Usually, $R^{2}>0.5$ reflects a significant relationship between $\boldsymbol{x}$ and $\boldsymbol{y}$.

4 Case study

4.1 Problem setup

To demonstrate the effectiveness and efficiency of the proposed surrogate, we consider the following benchmark Mindlin–Reissner model on the unit square domain Mindlin (\APACyear1951); Reissner (\APACyear1945):

[TABLE]

where $\Omega=[0,1]^{2}\subset\mathbb{R}^{2}$ is a smooth domain with a small thickness $t=0.1$, $\mathbf{C}$ is the positive definite tensor denoting the bending moduli, $\varepsilon$ is the linear Green strain tensor, $\boldsymbol{\theta}=[\theta_{x},\theta_{y}]$ represents the rotations of the surface, $w$ is the transverse displacement in z-direction, $\gamma$ denotes the scaled shear stresses, $f$ is the applied scaled transversal load, and $\lambda=E\kappa/2(1+v)$ is the shear modulus with $E$ denoting the Young’s modulus, $v=0.3$ denoting the Poisson ratio, and $\kappa=5/6$ denoting the shear correction factor. Essentially, Eq. 22 is a system of second order partial differential equations that can describe the bending behavior of a clamped plate loaded by a transverse force and Fig. 4. (a) shows the computational domain of our interest.

The objective of this benchmark problem is to build a surrogate model that can propagate and quantify the spatially-varying uncertainty associated with model inputs. The selected random input property is modeled by means of random field theory where its logarithm is a Gaussian random field:

[TABLE]

where $m(s)$ and $k\left(s,s^{\prime}\right)$ are the mean and covariance functions of the Gaussian process. In this case study, the mean function is defined to be zero and the exponentiated quadratic covariance function that is also known as the RBF kernel is adopted:

[TABLE]

where the correlation length $l=0.5$. We expand the scope of the case study by applying the proposed approach to the modeling of different mapping relationships, even though there is no obvious structure sharing between the input and output fields. Quantities of interest in different cases are summarized in Table 1:

4.2 Training, testing, and validation datasets

To obtain input-output data to construct the deep convolutional neural networks, we solve the forward PDEs using the finite element method (FEM) Hughes (\APACyear2012). After the convergence study (See Fig. 4. (b)), the unit square domain is discretized by $64\times 64$ quadrilateral isoparametric element (Q4) elements, ensuring the accuracy of the FE solution. Hence, the number of degrees of freedom of the high-fidelity FE model is $12675$. The FEM writes the weak form of equilibrium equations (Eq. 22) of the system into a set of algebraic equations, and solves the whole set by the Newton-Raphson method. The input to the FEM solver is the discretized representation of the random field(fields). For instance, each sample $\boldsymbol{x}_{i}\in\mathbb{R}^{64\times 64\times 2}$ in the case $3$ as the input covers the material and the load. The last column of Table 1 shows the dimensionality of the input-output relationship of each case.

There are several sampling schemes available to generate data that includes three sets: the training set, testing set, and validation set. We use the Latin hypercube sampling method to generate samples from the prescribed probability distributions. Specifically, we generate $1024$ input samples at which the high-fidelity model is evaluated, yielding the training dataset $\mathcal{D}_{train}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{1024}$. And we randomly sample another $200$ inputs for the testing $\mathcal{D}_{test}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{200}$. The reference statistics of the uncertainty analysis is computed by the Monte Carlo method using $10^{5}$ samples. Once the constructed network is optimally trained, we carry out the uncertainty quantification on the easy-to-evaluate surrogate instead of the computationally intensive finite element model with the same input snapshots.

The surrogate modeling algorithm is implemented in the open source machine learning library TensorFlow and the scripts are tested on a single NVIDIA GeForce GTX 1080 Ti X GPU. For the purpose of reproducibility, the code to replicate the optimization results in this work will be made available at https://xihaier.github.io upon publication of this manuscript.

4.3 Results

4.3.1 Case 1: one-to-one mapping

The first case concerns the use of field regressor to learn one-to-one mapping. Let the input be the material field, i.e. the young modulus and the output be the displacement field. The model inputs $\boldsymbol{x}\in\mathbb{R}^{4096}$ are discretized random field realizations and the model outputs $\boldsymbol{y}\in\mathbb{R}^{4096}$ are the transverse deformation at the finite element centers. The uncertainty analysis of the target system with high-dimensional stochastic inputs covers two parts. First, we have to train a surrogate model $\hat{f}:\mathbb{R}^{4096}\times[0,1]^{2}\rightarrow\mathbb{R}^{4096}$, which can accurately map the snapshot of the material field to the displacement field governed by the PDEs. Next, the trained surrogate is applied to the uncertainty quantification tasks, computing statistical moments and PDFs.

Network architecture. After extensive hyperparameter and network architecture search, the field regressor is built using 21 convolutional layers and this configuration is hence referred to as FR-21. Specifically, the dimension of each random field sample has been reduced twice by the downsampling processor $64\times 64\rightarrow 31\times 31\rightarrow 17\times 17$. There are three dense blocks within FR-21. They share the same kernel size, stride number, and padding number while have different growth number and depth. Fig. 5 gives a general summary of the one-to-one mapping problem, graphic illustration of the data flow of FR-21, and the details of these network configuration parameters.

To learn the high-dimensional input-output relationship, the gradient-based optimizer Adam is selected. We train the surrogate model (FR-21) with training data size varies from $64$ to $1024$. The initial learning rate is set to $\eta_{0}=0.005$ and the step decay learning strategy (often known as annealing) is adopted to gradually lower the learning rate during training. The default value $\eta_{0}$ is reduced by a annealing rate $\zeta=0.75$ every $20$ epochs and the model trained with $500$ epochs in total. For regularization, the weight decay is initialized to $\lambda=7\times 10^{-6}$ and the batch size is $8$.

Optimization results. Fig. 6 (A) plots the training process in terms of the RMSE value, which is recorded every $20$ epochs. The rapidness and steadiness of the convergence of the FR-21 can be seen from the results. In particular, after around $300$ epochs of the Adam optimizer, an approximately stabilized solution is achieved. The optimization algorithm converges faster with more training samples and the RMSE value of the testing dataset reduces as the training size increases: $64\rightarrow\dots\rightarrow 1024$. Meanwhile, the coefficient of determination (also known as the $R^{2}$-score) is evaluated using 200 test samples with training dataset sizes. Note that the trained FR-21 shows a remarkable generalization ability where the $R^{2}$-score is more than $0.9$ in all the trials. In Fig. 6 (B), the model achieves a relatively high $R^{2}$-score of $0.967$ with only $64$ training samples and the $R^{2}$-score gets close to $0.99$ when increasing the training dataset to $512$ samples, which indicates that the predicted output fields by FR-21 match the ground truth computed by the FEM very well. Fig. 6 (C) show some of the predicted results. The errors of the predicted solutions of FR-21 are compared with respect to the FE solutions. Note that the results representing the output fields have been standardized so that a more clear comparison can be carried out between FR-21 with different training data sizes. For visualization purpose, we randomly selected a sample from the testing dataset $\mathcal{D}_{test}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{200}$. Remarkably, the propsoed FR-21 is capable of providing good predictions. Even when using only $64$ training samples, the model is able to quickly capture the landscape of the target output field. The prediction quality can be significantly improved as more samples are put into the training process.

Uncertainty analysis results. The uncertainty analysis results predicted by the FR-21 is verified against the solution computed by the FEM solver. The number of training samples is set to $N=1024$, and the number of testing samples is set to $M=1\times 10^{5}$. Fig. 6 (A) and (B) show the the mean and variance of the displacement field, respectively. It can be seen that FR-21 achieves smooth contour estimation that is close to the Monte Carlo result. The maximum error of the estimate measured in a relative error form is less than $10\%$. Meanwhile, we randomly select two points out of the domain and show the PDFs, respectively. The unbounded kernel smoothing function is used to estimate the PDFs, where the type of kernel smoother is set to normal. Fig. 6 (C) compares the PDFs at two locations from FR-21 and FEM solver. As desired, the PDFs obtained by the surrogate model at two randomly selected locations are almost identical to the reference solution.

4.3.2 Case 2: one-to-many mapping

In this case, the field regressor is generalized to the surrogate learning of one-to-many mapping. The input field is the young modulus $\boldsymbol{x}=E(s)$, and the output fields are stress fields that essentially are functions of displacement fields $\boldsymbol{y}=\left[\sigma_{v}(s),\tau_{\max}(s),\tau_{xy}(s)\right]$ (See Table 1 for details). The objective is to build a surrogate model that can predict three discretized random fields simultaneously from one input random field $\hat{f}:\mathbb{R}^{64\times 64}\times[0,1]^{2}\rightarrow\mathbb{R}^{64\times 64\times 3}$.

Network architecture. The surrogate model is derived using 25 convolutional layers in total, and hence is referred to as FR-25. There are four dense blocks in the FR-25. In particular, we build a dense block before the last convolutional layer, ensuring that sufficient information has been extracted to construct three independent output fields. The number of feature maps increases from $72$ before the block to $144$ after a set of block layers $21-24$. Fig. 8 gives a graphic representation of the one-to-many mapping problem $\mathbb{R}^{4096}\rightarrow\mathbb{R}^{12288}$ and summarizes the detailed designs of FR-25.

Optimization results. The surrogate model is trained with the same optimization configuration stated in Section 4.3.1. Fig. 9 (A) confirms the learning effectiveness of our proposed surrogate model. The Adam optimizer reaches stabilized state after $350$ epochs. The objective function minimizes faster with more training data, and the RMSE value at epoch $500$ increases when the training size decreases. Fig. 9 (B) shows the $R^{2}$-score for each training trial. Note that the $R^{2}$-score is evaluated using all the predicted fields, implying each sample $\hat{\boldsymbol{y}}_{i}$ in Eq. 21 belongs to the real-valued vector space $\mathbb{R}^{64\times 64\times 3}$. The results show that the surrogate, FR-25, provides an accurate approximation of the high-dimensional input-output relationship even with a few samples (e.g. 64 samples). Specifically, these $R^{2}$-score values indicate the majority of the observed variation (approximately 90 percents) can be explained by the model. We present the predicted results using $64$ and $1024$ training samples in Fig. 9 (C), respectively. The estimated stress fields via the FR-25 closely resemble the landscape of the ground truth data when only $64$ samples are available for the model training. And the prediction accuracy improves substantially when the training sample size increases to $1024$. Noticebly, the preiction resutls of the shear stress field $\tau_{xy}(s)$ is relatively less smoother, especially when the training sample size is small. This is due to the weights assigned to each output field. Without any prior knowledge, the intuitive assignment is the use of uniform weights. The different performances regarding different output fields are largely alleviated as we increase the training samples, and the predictions become very close to the reference solutions, where the relative maximum error of $\sigma_{v}(s)$, $\tau_{max}(s)$, and $\tau_{xy}(s)$ are in the same range.

Uncertainty analysis results. The predictive performance of the FR-25 to the machine learning with a small dataset seen in this case is a strong indication of the effectiveness and robustness of using the proposed surrogate modeling approach for uncertainty analysis. A number of $1\times 10^{5}$ independent and identically distributed samples have been generated for the computation of the reference statistics using the Monte Carlo method. With respect to the PDFs, we once again randomly select $2$ points from each output field and use kernel methods to estimate the corresponding PDF. For the mean field prediction, the maximum error is less than $5\%$ measured in a relative error sense and the error centers around the extreme value region, for instance, the central cross region of the $\tau_{max}(s)$ field (See Fig. 10 (A)). For the mean variance prediction, the maximum error is distributed on the boundary lines in terms of the prediction of $\sigma_{v}(s)$ and $\tau_{max}(s)$ and is overlapped with the extreme value region in the case of $\tau_{xy}(s)$ (See Fig. 10 (B)). The PDFs obtained by the surrogate model at two randomly selected locations are almost identical to the reference solution. Fig. 10 (C) shows the randomly selected points have non-gaussian distributions of the $\sigma_{v}(s)$, $\tau_{max}(s)$, and $\tau_{xy}(s)$ value. And the FR-25 can capture different PDFs with remarkably high precision.

4.3.3 Case 3: many-to-many mapping

The third and final case study investigates the efficacy of the proposed surrogate modeling approach for many-to-many mapping problems. We consider the boundary value problem defined by the elliptic PDEs (Eq. 22) with multiple random input fields. In particular, the material $E(s)$ and loading $f(s)$ properties are considered to vary spatially where the correlated variation has the same analytical form (Eq. 23). Subsequently, we are interested in characterizing the spatial variability induced effects on the structural performances that are denoted by the displacement $w(s)$ and stress $\sigma_{v}(s)$. The input-output relationship satisfies the high-dimensional mapping $\hat{f}:\mathbb{R}^{64\times 64\times 2}\times[0,1]^{2}\rightarrow\mathbb{R}^{64\times 64\times 2}$ on the unit square domain.

Network architecture. Similar to the previous case, the deep neural networks based surrogate model is configured with 25 convolutional layers. Since the model input $\boldsymbol{x}$ covers two independent random fields that are physically irrelevant, more features have been extracted and stacked by the end of the first layer, where the number of input and output features is $2$ and $64$, respectively. As displayed in the table of Fig. 11, a total of $4$ dense blocks have been constructed on account for the sufficiency of the reconstruction of the output fields. The concept of the many-to-many mapping problem is illustrated in Fig. 11 (A) and a graphic diagram of the FR-25 in the present case is given in Fig. 11 (B).

Optimization results. Fig. 12 (A) shows the optimization history of the model learning with different size of training samples. The strength of the proposed surrogate in terms of training is verified as the objective function decreases dramatically within the first $200$ epochs. In comparing the different trials of this case study, again, it can be seen that an increase in the training sample size leads to a better-optimized result. Indeed, the $R^{2}$-score illustrates the importance of providing sufficient training samples, where the coefficient of determination value is relatively lower compared to the second case with the same amount of training samples, with the $R^{2}$-score of the trial using $64$ samples slightly smaller than $0.9$. This is due to the fact that the complexity of this case is higher than the second one. Potential improvement strategy includes altering the architecture design of the field regressor, for instant, extracting initial features in a separate manner rather than the joined style, that is,

[TABLE]

For the purpose of easy comparison, we attempted to use a united design in all cases. The accuracy and consistency of the predicted results of the proposed approach can once again be seen from Fig. 12 (C), (D), and (E). In particular, the precision observed in the selected trials of this case study approximately equivalent to those seen in the previous cases, and thus indicates the robustness and strong potential of the proposed surrogate modeling approach for solving a wide variety of problems of interest.

Uncertainty analysis results. Fig. 13 first investigates the prediction performance of the mean and variance fields by means of FR-25. The reference solution is obtained by numerically solving the governing SPDEs with the finite element scheme $10^{5}$ times. Fig. 13 (A) and (B) present the good agreement between the surrogate modeling and high-fidelity FE results for the statistics of the target output fields. Fig. 13 (C) also depicts the PDFs of two output fields at three specific locations. It can be seen that our field regressor is able to accurately approximate the probabilistic response. This once again illustrates the robustness and efficiency of the proposed field regressor, i.e., using $1024$ samples is sufficient to capture most of the stochasticity in this $\mathbb{R}^{8192}\rightarrow\mathbb{R}^{8192}$ SPDEs.

5 Conclusion

This paper presents a deep neural network approach for uncertainty analysis of elliptic systems with spatially varying properties. Note that conventional random variable representation with uniform distribution of variance cannot consider the effect of spatial randomness notwithstanding the fact that spatial variability is naturally associated with a wide range of engineering systems. In the present study, we describe the correlation that exists in the uncertain spatial distribution by means of random field theory. After spatial discretization, the stochastic parameter space in terms of the modeling of the input-output relationship is high-dimensional. To overcome this high-dimensional stochasticity issue, we introduce a hierarchical surrogate model named the field regressor utilizing convolutional neural networks as the model basis. Compared to conventional approaches, our field regressor provides a direct approximation of the intrinsic input/output relationships without mapping the high-dimensional data to a lower dimension.

Specifically, novelties of the proposed field regressor can be outlined in three perspectives. First, information flowing over the model is gradually downsized to a stack of feature maps. Model-estimated output fields are then reconstructed from extracted features. Such a multi-scaling data mechanism hereby advances the training process as fewer model parameters are required. Secondly, short connections between nonadjacent convolutional layers are introduced to the basic network architecture, allowing previously extracted features to be reused in the subsequent layers. This property favors the machine learning process as a deeply structured architecture becomes feasibly trainable and numerically stable. Thirdly, we use bicubic resizing and convolution operation in the upsampling process. Consequently, detailed variation patterns in the high-dimensional data can be better preserved.

The effectiveness and efficiency of the proposed field regressor is demonstrated on a benchmark UQ problem. It is observed that the field regressor is capable of directly inferring a variety of high-dimensional mapping relationships including one-to-one ($\mathbb{R}^{4096}\rightarrow\mathbb{R}^{4096}$), one-to-many ($\mathbb{R}^{4096}\rightarrow\mathbb{R}^{12288}$), and many-to-many ($\mathbb{R}^{8192}\rightarrow\mathbb{R}^{8192}$) mappings. Even when there is no obvious functional relationship between the data structures of the input and output fields, the proposed field regressor can accurately capture the intrinsic mapping relationships. The use of the field regressor greatly accelerates the uncertainty analysis, i.e. propagating the random field represented uncertainty to the output quantities of interest.

It should be noted that the current surrogate is merely built through data without integrating any physics into the neural network architecture. Possible improvement strategies of the model performance may include utilizing well-established physical interpretations during the architecture design process of the surrogate. For instance, building a new neural network as a representation of governing equations or engineering constraints in addition to the current field regressor will be helpful. In the future, besides applying the field regressor to other potential engineering applications with high-dimensional data, the current surrogate modeling technique will be extended to spatiotemporal problems.

Acknowledgements.

X.L. acknowledges the NatHaz Modeling Laboratory at the University of Notre Dame where this work was performed. This research effort was supported in part by the National Science Foundation through Grant No. 1520817 and No. 1612843. The codes and data used in this work will be made available at https://xihaier.github.io/ upon publication of this manuscript.