Bayesian deep learning with hierarchical prior: Predictions from limited and noisy data

TL;DR

This paper introduces a Bayesian deep learning framework with hierarchical priors that effectively handles limited and noisy data, providing robust predictions and uncertainty quantification through variational inference and Monte Carlo methods.

Contribution

It integrates hierarchical Bayesian modeling with deep learning, enabling automatic uncertainty quantification and improved robustness in noisy, limited data scenarios.

Findings

Enhanced predictive accuracy with contaminated data

Effective uncertainty quantification via Monte Carlo simulation

Robustness of hierarchical Bayesian approach demonstrated

Abstract

Datasets in engineering applications are often limited and contaminated, mainly due to unavoidable measurement noise and signal distortion. Thus, using conventional data-driven approaches to build a reliable discriminative model, and further applying this identified surrogate to uncertainty analysis remains to be very challenging. A deep learning approach is presented to provide predictions based on limited and noisy data. To address noise perturbation, the Bayesian learning method that naturally facilitates an automatic updating mechanism is considered to quantify and propagate model uncertainties into predictive quantities. Specifically, hierarchical Bayesian modeling (HBM) is first adopted to describe model uncertainties, which allows the prior assumption to be less subjective, while also makes the proposed surrogate more robust. Next, the Bayesian inference is seamlessly integrated into the DL framework, which in turn supports probabilistic programming by yielding a probability distribution of the quantities of interest rather than their point estimates. Variational inference (VI) is implemented for the posterior distribution analysis where the intractable marginalization of the likelihood function over parameter space is framed in an optimization format, and stochastic gradient descent method is applied to solve this optimization problem. Finally, Monte Carlo simulation is used to obtain an unbiased estimator in the predictive phase of Bayesian inference, where the proposed Bayesian deep learning (BDL) scheme is able to offer confidence bounds for the output estimation by analyzing propagated uncertainties. The effectiveness of Bayesian shrinkage is demonstrated in improving predictive performance using contaminated data, and various examples are provided to illustrate concepts, methodologies, and algorithms of this proposed BDL modeling technique.

Tables

Table 1. Table 1: Comparison of the coefficient of determination ( R 2 superscript 𝑅 2 R^{2} ) of the different surrogate models where the training dataset is contaminated by different noise levels.

Method	clean data	$𝝈=0.1$	$𝝈=0.3$	$𝝈=0.5$	$𝝈=0.7$	$𝝈=0.9$
PR	0.9950	0.9937	0.9884	0.9807	0.9680	0.9533
SVM	0.9783	0.9784	0.9758	0.9604	0.9486	0.9388
NN	0.9890	0.9854	0.9828	0.9770	0.9526	0.9497
BDL	0.9883	0.9893	0.9928	0.9757	0.9672	0.9516

Table 2. Table 2: Comparison of the root mean squared error (RMSE) of the different surrogate models where the training dataset is contaminated by different noise levels.

Method	clean data	$𝝈=0.1$	$𝝈=0.3$	$𝝈=0.5$	$𝝈=0.7$	$𝝈=0.9$
PR	0.2483	0.2629	0.3917	0.4863	0.6227	0.7643
SVM	0.5397	0.5424	0.5642	0.6997	0.7934	0.9198
NN	0.3817	0.4068	0.4783	0.5276	0.7350	0.8077
BDL	0.3270	0.3098	0.2964	0.5425	0.6091	0.7933

Table 3. Table 3: Statistics of the uncertain input parameters for the flat membrane.

Basic variables	First parameter	Second parameter	Distribution type
$l$	$\mu =2$	$\sigma =0.05$	Normal
$b$	$\mu =1$	$\sigma =0.05$	Normal
$t$	$min=0.001$	$max=0.002$	Uniform
$E$	$\mu =210$	$\sigma =10$	Normal
$\upsilon$	$\mu =\log \left(0.3\right)$	$\sigma =0.01$	Lognormal

Table 4. Table 4: Comparison of the coefficient of determination.

Prior type	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{64}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{128}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{256}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{512}$
Direct prior	0.9720	0.9982	0.9990	0.9993
Hyper prior	0.9935	0.9989	0.9983	0.9994

Table 5. Table 5: Comparison of the root mean square error.

Scenario	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{100}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{200}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{300}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{400}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{500}$	${\{{𝒙}_i,{𝒚}_i\}}_{i=1}^{600}$
${\text{Windward}}^{†}$	0.077	0.030	0.027	0.022	0.019	0.017
${\text{Windward}}^{‡}$	0.051	0.021	0.021	0.020	0.018	0.015
${\text{Leftside}}^{†}$	0.039	0.030	0.027	0.019	0.013	0.012
${\text{Leftside}}^{‡}$	0.037	0.030	0.024	0.016	0.012	0.012
${\text{Leeward}}^{†}$	0.014	0.013	0.011	0.011	0.010	0.010
${\text{Leeward}}^{‡}$	0.014	0.011	0.011	0.010	0.010	0.009
${\text{Rightside}}^{†}$	0.037	0.025	0.025	0.021	0.019	0.015
${\text{Rightside}}^{‡}$	0.033	0.025	0.023	0.020	0.014	0.011

Equations

$$\mathbf{y}=\mathcal{F}\left(\mathbf{x}\right)\stackrel{\approx }{\to }\hat{\mathbf{y}}=\hat{\mathcal{F}}\left(\mathbf{x}\right)$$

$$\hat{\mathcal{F}}\left(\mathbf{x}\right)={\hat{f}}^K\circ {\hat{f}}^{K-1}\circ \dots {\hat{f}}^1\left(\mathbf{x}\right)$$

$${\hat{f}}^i\circ {\hat{f}}^j\doteq {\hat{f}}^i\left({\hat{f}}^j\left(\cdot \right)\right)$$

$${\mathbf{z}}^i={\hat{f}}_1^i\left({\mathbf{x}}^i\right)={\mathbf{\omega }}^i{\mathbf{x}}^i+{\mathbf{b}}^i$$

$${\hat{\mathbf{y}}}^i={\hat{f}}_2^i\left({\mathbf{x}}^i\right)=\sigma \left({\mathbf{z}}^i\right)$$

$${\mathbf{\omega }}^{\star }=\mathrm{argmin}\sum _{i=1}^n\mathcal{L}\left({\mathbf{y}}_i,\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)$$

$$\mathcal{L}\left({\mathbf{y}}_i,\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)=p\left(\mathcal{D}|\mathbf{\omega }\right)=p\left(\mathbf{y}|\mathbf{x},\mathbf{\omega }\right)$$

$$\mathcal{L}\left({\mathbf{y}}_i,\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)$$

$$={\left(\frac{1}{{\left(2\pi {\mathbf{\sigma }}_{\mathbf{ϵ}}^2\right)}^n}\right)}^{\frac{1}{2}}\prod _{i=1}^n\exp \left(-\frac{1}{2{\mathbf{\sigma }}_{\mathbf{ϵ}}^2}{\left({\mathbf{y}}_i-\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)}^2\right)$$

$${\mathbf{\omega }}^{\star }=\mathrm{argmin}\sum _{i=1}^n\left(-\frac{1}{2{\mathbf{\sigma }}_{\mathbf{ϵ}}^2}{\left({\mathbf{y}}_i-\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)}^2\right)-\frac{n}{2}\log \left(2\pi {\mathbf{\sigma }}_{\mathbf{ϵ}}^2\right)$$

$${\tau }_{\mathbf{ϵ}}^{MLE}={\left(\frac{1}{n{N}_{\mathbf{\omega }}}\sum _{i=1}^n{\left({\mathbf{y}}_i-\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)}^2\right)}^{-1}$$

$${\mathbf{\omega }}^{\star }=\mathrm{argmin}\sum _{i=1}^n{\left({\mathbf{y}}_i-\hat{\mathcal{F}}\left({\mathbf{x}}_i\right)\right)}^2$$

$${\mathbf{\omega }}_{t+1}\leftarrow {\mathbf{\omega }}_t+\eta \nabla \mathcal{L}\left(\mathbf{\omega }\right)$$

$$\nabla \mathcal{L}\left(\cdot \right)=\mathbb{E}[g\left(\mathbf{x}\right)]$$

$$g\left(\mathbf{x}\right)=\frac{{\hat{M}}_t}{\sqrt{{\hat{V}}_t}+ϵ}=\frac{M_t}{1-{\beta }_1^t}/\left(\sqrt{\frac{V_t}{1-{\beta }_2^t}}+ϵ\right)$$

$$M_t$$

$$V_t$$

$$\mathbf{\omega }\sim p\left(\mathbf{\omega }\right)$$

$$p\left(\mathbf{\omega }\right)=\int p\left(\mathbf{\omega },\mathbf{\eta }\right)d\mathbf{\eta }$$

$$p\left(\mathbf{\omega },\mathbf{\eta }\right)=p\left(\mathbf{\omega }|\mathbf{\eta }\right)p\left(\mathbf{\eta }\right)$$

$$p\left(\mathbf{\omega }\right)=\prod _{i=1}^{K}p\left({\mathbf{\omega }}_i|{\mathbf{\mu }}_{{\mathbf{\omega }}_i},{\tau }_{{\mathbf{\omega }}_i}\right)=\prod _{i=1}^K\mathcal{N}\left({\mathbf{\mu }}_{{\mathbf{\omega }}_i},{\tau }_{{\mathbf{\omega }}_i}^{-1}\mathbf{I}\right)$$

$$p\left(\mathbf{\omega }|{\tau }_{\mathbf{\omega }},{\alpha }_{\tau },{\beta }_{\tau }\right)=p\left(\mathbf{\omega }|{\tau }_{\mathbf{\omega }}\right)p\left({\tau }_{\mathbf{\omega }}|{\alpha }_{\tau },{\beta }_{\tau }\right)\text{where}{\alpha }_{\tau }>0,{\beta }_{\tau }>0$$

$$p\left(\mathbf{\omega }|{\alpha }_{\tau },{\beta }_{\tau }\right)$$

$$p\left(\mathbf{ϵ}|{\mathbf{\mu }}_{\mathbf{ϵ}},{\mathbf{\sigma }}_{\mathbf{ϵ}}^2\right)=\mathcal{N}\left({\mathbf{\mu }}_{\mathbf{ϵ}},{\mathbf{\sigma }}_{\mathbf{ϵ}}^2\right)$$

$$\mathbf{y}\sim \mathcal{N}\left(\mathbf{y}|{\mathbb{E}}_{p\left(\mathbf{\omega }|\mathcal{D}\right)}[\hat{\mathcal{F}}\left(\mathbf{x}\right)],{\tau }_{\mathbf{ϵ}}\mathbf{I}\right)$$

$$p\left(\mathbf{\omega },\mathbf{\eta }|\mathcal{D}\right)=\frac{p\left(\mathcal{D}|\mathbf{\omega },\mathbf{\eta }\right)p\left(\mathbf{\omega }|\mathbf{\eta }\right)p\left(\mathbf{\eta }\right)}{p\left(\mathcal{D}\right)}$$

$$p\left(\mathcal{D}\right)={\int }_{\mathbf{\omega }}p\left(\mathcal{D}|\mathbf{\omega }\right)p\left(\mathbf{\omega }\right)d\mathbf{\omega }$$

$$p\left(\mathbf{\omega }|\mathcal{D}\right)\approx q_{\mathbf{\xi }}\left(\mathbf{\omega }\right)\in \mathbf{\Xi }$$

$$p\left(\mathbf{\omega }|\mathcal{D}\right)\simeq q_{{\mathbf{\xi }}^{\star }}\left(\mathbf{\omega }\right)$$

$$=\mathrm{argmin}\int q_{\mathbf{\xi }}\left(\mathbf{\omega }\right)\log \frac{q_{\mathbf{\xi }}\left(\mathbf{\omega }\right)}{p\left(\mathbf{\omega }|\mathcal{D}\right)}d\mathbf{\omega }$$

$$q_{{\mathbf{\xi }}^{\star }}\left(\mathbf{\omega }\right)$$

$$=\mathrm{argmax}{\mathbb{E}}_{q_{\mathbf{\xi }}\left(\mathbf{\omega }\right)}\left[\log p\left(\mathcal{D}|\mathbf{\omega }\right)\right]-\mathbb{KL}\left(q_{\mathbf{\xi }}\left(\mathbf{\omega }\right)|p\left(\mathbf{\omega }\right)\right)$$

Full text

††journal: Elsevier

1 Introduction

Applications of data-driven approaches for learning the performance of engineering systems using limited experimental data is hindered by at least three factors. First, the original input-output patterns are often governed by a series of highly nonlinear and implicit partial differential equations (PDEs), and hence approximation of their functional relationships may be proportionally computation demanding [1, 2]. Secondly, interpolation and extrapolation techniques are usually needed to extract knowledge from acquired data in consideration of only a limited number of sensors used in practice along with the fact that sensor malfunction often occurs in real time. Nevertheless, it is very difficult to establish an accurate discriminative model merely from data, especially when a relatively small dataset is available [3, 4]. Thirdly, experimental data is inevitably contaminated by noise from different sources, for example, signal perturbation induced noisy sensing during monitoring. The performance of conventional discriminative algorithms may be noticeably impaired if proper noise reduction has not been performed [5, 6]. In this context, we present a machine learning based predictive model that is capable of providing high-quality predictions from limited and noisy data.

To date, most machine learning models are deterministic, which implies that a certain input sample $\boldsymbol{x}_{i}$ is strictly bounded to a point estimator $\hat{\boldsymbol{y}}_{i}$ notwithstanding the existence of model uncertainties [7, 8]. Probabilistic modeling, on the other hand, emerges as an attractive alternative on account of its ability to quantify the uncertainty in model predictions, which can prevent a poorly trained model from being overconfident in predictions, and hence helps stakeholders make a more reliable decision [7, 9]. In literature, Gaussian processes (GPs) and generalized polynomial chaos (gPC) are two representative members from the probabilistic modeling family [9, 10, 11, 12]. From a mathematical standpoint, GPs put a joint Gaussian distribution over input random variables by defining a mean function $E[\boldsymbol{x}]$ and a covariance function $Cov[\boldsymbol{x}_{i},\boldsymbol{x}_{j}]$ [10]. Then, GPs compute hyperparameters of the spatial covariance function and propagte the inherent randomness of $\boldsymbol{x}$ in virtue of the Bayes’ theorem. Next, gPC is an effective way to propagate uncertain quantities by means of utilizing a set of random coefficients $\{\beta_{1},\beta_{2},\dots,\beta_{n}\}$ and orthogonal polynomials $\{\phi_{1},\phi_{2},\dots,\phi_{n}\}$. Approximation approaches (e.g. Galerkin projection) are usually used to determine the unknown coefficients $\boldsymbol{\beta}$ and polynomial basis $\boldsymbol{\phi\left(\boldsymbol{x}\right)}$ [12]. Even though GPs and gPC are capable of computing empirical confidence intervals, the inference complexity may become overwhelming when the number of observations increases (e.g. cubic scaling relationship $\mathcal{O}\left(N^{3}\right)$ between the computational complexity and data $\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{N}$ is found in the GPs case [13]). Furthermore, scaling a GPs model or identifying random coefficients $\boldsymbol{\beta}$ of a gPC model for problems with high-dimensional data remains challenging [3, 9, 10].

On the contrary, deep learning (DL) has differentiated itself within the realm of machine learning for its superior performance in handling large-scale complex systems. With the strong support of high-performance computing (HPC), DL has made significant accomplishments in a wide range of applications such as image recognition, data compression, computer vision, and language processing [14, 15]. Nonetheless, the same level of application has not been observed for its probabilistic version [16, 17, 18]. This paper attempts to bridge the modeling gap between DL and Bayesian learning by presenting a new paradigm Bayesian deep learning (BDL) model. With the aim of developing a surrogate model that can be used to accelerate the uncertainty analysis of engineering systems using noisy data, we focus on three main aspects in probabilistic modeling (See Fig. 1).

First, Bayesian statistical inference encodes the subjective beliefs, whereas prior distributions are imposed on model parameters to represent the initial uncertainty. A conventional DL supported surrogate $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$, however, can be captivating and confusing in equal measure at the same time. Owning a deeply structured network architecture allows $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ to approximate a wide variety of functions, but also makes model parameters hard to interpret, which in turn increases the difficulty of choosing a reasonable prior distribution for $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ [19, 20]. In [19], Lee shows noninformative prior (e.g. Jeffreys prior) that bears objective Bayes properties is liable to be misled by the variability in data. And using the Fisher information matrix to compute a Jeffreys prior for a large network architecture can be computationally prohibitive [21]. On the informative prior side, zero-centered Gaussian is, not surprisingly, extensively explored in early work on Bayesian neural networks (BDL) on account of its flexibility in implementation as well as its natural regularization mechanism [16, 17], where the quadratic penalty for BDL parameters alleviates the overfitting problem. Later in [18], Neal points out that employing a heavy-tailed distribution (e.g. Cauchy distribution) to represent prior knowledge can provide a more robust shrinkage estimator and diminish the effects of outlying data. Nonetheless, Cauchy distribution is difficult to implement because it does not have finite moments [20]. To develop a prior model that is more amenable to reform and work with, we investigate the efficacy of applying hierarchical Bayesian modeling (HBM) to define the prior distribution. And it is found that HBM can efficiently ameliorate the prior assumptions induced model performance variance by diffusing the influences that are cascaded down from the top level, hence allowing a relatively more robust distribution model.

Secondly, the determination of the posterior distribution $p\left(\boldsymbol{\omega}|\boldsymbol{\mathcal{D}}\right)$ requires integrating model parameters $\{\omega_{1},\omega_{2},\dots,\omega_{n}\}$ out of the likelihood function. Unfortunately, performing numerical integration in BDL’s parameter space is always computationally intractable as a neural networks model is commonly configured with hundreds/thousands of parameters [22, 23]. Approximation methods that can be broadly classified into the sampling method and the optimization method have been introduced to alleviate such computational bottleneck [24]. For the the sampling method, Markov Chain Monte Carlo (MCMC) has been explored in early work to calculate the posterior probabilities [17, 20, 25]. The main idea of MCMC is to produce numerical samples from the posterior distribution by simulating a discrete but dependent Markov chain $\{\boldsymbol{\omega}_{i}\}_{i=1}^{M}$ on the state space, where $\pi\left(\boldsymbol{\omega}\right)\approx p\left(\boldsymbol{\omega}|\boldsymbol{\mathcal{D}}\right)$. The sufficient condition that ensures the stationary distribution converges to the target posterior distribution requires the transition kernel $T\left(\cdot\right)$ to have detailed balance properties $p\left(\boldsymbol{\omega}|\boldsymbol{\mathcal{D}}\right)T\left(\boldsymbol{\omega}^{{}^{\prime}}\rightarrow\boldsymbol{\omega}\right)=p\left(\boldsymbol{\omega}^{{}^{\prime}}|\boldsymbol{\mathcal{D}}\right)T\left(\boldsymbol{\omega}\rightarrow\boldsymbol{\omega}^{{}^{\prime}}\right)$. However, the convergence of MCMC algorithms can be extremely slow in the presence of large datasets since the burn-in process to eliminate the initialization bias is greatly extended [9, 23]. More recently, variational inference (VI) method has been employed for inferring an intractable posterior distribution where the probabilistic inference problem is cast in a deterministic optimization form [26, 27, 28]. A proxy probability distribution $q\left(\boldsymbol{\omega}\right)$, which is formally explicit and computationally efficient, is introduced to approximate the true posterior distribution $p\left(\boldsymbol{\omega}|\boldsymbol{\mathcal{D}}\right)$. Compared to the sampling method, VI has the advantage of approximating non-conjugate distributions by virtue of optimizing an explicit objective function [23, 24], and VI solves the intractable integrals in a more efficient manner on account of off-the-peg optimization algorithms can be seamlessly adapted to the minimization problem [26, 28, 29]. Unfortunately, the differentiation of the objective function regarding the proxy posterior involves the determination of the expectations with respect to the variational parameters where Monte Carlo gradient estimator may give a high variance [30, 31]. To address this issue, we reparameterize our objective function by introducing a set of auxiliary variables. It should be noted that such reparameterization would not only yield an unbiased estimator of the objective function but also provides an efficient approximation of variational parameters via permitting the use of stochastic gradient descent (SGD) in optimization.

Lastly, Monte Carlo (MC) method is used in the predictive phase of Bayesian inference. MC method draws numerical samples from the proxy probability distribution $q\left(\boldsymbol{\omega}\right)$, builds a predictive probability distribution of new data, and assigns a confidence level to the model prediction for representing model uncertainty. The following outline of this paper is intended as: Section 2 gives a brief introduction of surrogate modeling using deep neural networks. Section 3 describes the proposed BDL model in detail. In Section 4, various examples are provided to demonstrate the effectiveness of BDL in dealing with noisy data. Finally, Section 5 draws major conclusions.

2 Deterministic modeling: a deep learning framework

2.1 Neural networks based surrogate model

In the context of supervised learning [9, 15], let $\boldsymbol{x}=[x_{1},x_{2},\dots,x_{m}]\in\mathbb{R}^{m}$ denote an input vector, and the corresponding output vector $\boldsymbol{y}=[y_{1},y_{2},\dots,y_{n}]\in\mathbb{R}^{n}$ is estimated by a computationally intensive model $\mathcal{M}\left(\cdot\right)$ (e.g. a large finite element model). We are interested in using neural networks to approximate functional relationships $\mathcal{F}\left(\cdot\right)$ between $\boldsymbol{x}$ and $\boldsymbol{y}$.

[TABLE]

where $\hat{\mathcal{F}}\left(\cdot\right)$ is the mathematical expression of neural networks based surrogate $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$, and theoretically it can be proportionately broken down as:

[TABLE]

with $K$ denoting the layer number, and $\circ$ symbolizing the functional composition operation which is defined as [32]:

[TABLE]

Each function in the sequence $\hat{f}^{i}\left(\cdot\right),i=1,2,\dots,K$ contains two steps, where the first step is identical to a linear regression:

[TABLE]

In Eq. 4, $\boldsymbol{x}^{i}$ is the input vector of the $i^{th}$ layer, $\boldsymbol{\omega}^{i}$ is the weight matrix, and $\boldsymbol{b}^{i}$ is the $i^{th}$ bias term. For the sake of brevity, $\boldsymbol{b}^{i}$ can be integrated into $\boldsymbol{\omega}^{i}$ by introducing an additional input variable $x_{0}^{i}=1$. For the rest of the paper, let $\boldsymbol{\omega}$ be a tensor containing all model parameters [32, 33]. Next, $\hat{f}^{i}\left(\cdot\right)$ applies an element-wise nonlinear transformation to the intermediate output $\boldsymbol{z}^{i}$ in the second step:

[TABLE]

where $\sigma\left(\cdot\right)$ is often referred to as the activation function [14, 15]. Selection of $\sigma\left(\cdot\right)$ directly depends on the characteristics of $\mathcal{M}\left(\cdot\right)$. Moreover, a network architecture is deemed to be deep when $K>3$ [14]. Thus, a deep learning framework can be effectively built by increasing the composition size $K$.

2.2 Probabilistic interpretation of $L^{2}$ loss function

Consider a parameterised deep neural network model $\hat{\boldsymbol{y}}=\hat{\mathcal{F}}_{\boldsymbol{\omega}}\left(\boldsymbol{x}\right)$ described in Section 2.1 and a training dataset $\boldsymbol{\mathcal{D}}=\{\boldsymbol{x}_{i},\boldsymbol{y}_{i}\}_{i=1}^{N}$, the next step is to find an optimal $\boldsymbol{\omega}^{\star}$ such that the surrogate $\hat{\mathcal{F}}_{\boldsymbol{\omega}}\left(\boldsymbol{x}\right)$ best describes the data $\boldsymbol{\mathcal{D}}$. In this regard, a loss function $\mathcal{L}\left(\cdot\right)$ that measures the error between the predicted value $\hat{\boldsymbol{y}}$ and the expected result $\boldsymbol{y}$ is defined. Then, $\hat{\mathcal{F}}_{\boldsymbol{\omega}}\left(\boldsymbol{x}\right)$ is trained to approximate $\mathcal{F}\left(\boldsymbol{x}\right)$ by minimizing the empirical loss through tuning model parameters $\boldsymbol{\omega}$:

[TABLE]

From a probabilistic modeling perspective, the interest of loss function $\mathcal{L}\left(\cdot\right)$ is in the probability distribution of $\boldsymbol{y}$ as a funciton of $\boldsymbol{x}$:

[TABLE]

where the model parameters can be learned by the method of maximum likelihood estimate (MLE) [9, 22], which searches an estimator for $\boldsymbol{\omega}$ that maximizes the likelihood term $\boldsymbol{\omega}^{MLE}=\operatorname*{arg\,max}p\left(\boldsymbol{\mathcal{D}}|\boldsymbol{\omega}\right)$. Let noise term $\boldsymbol{\epsilon}$ be independent and identically distributed (i.i.d.), the probability density associated with paired observations under Gaussian assumption can be expressed as:

[TABLE]

Usually, the numerical implementation of the MLE method performs the minimization problem of Eq. 6 in a logarithmic scale:

[TABLE]

where the precision term $\tau_{\boldsymbol{\epsilon}}$ is determined by minimizing the negative log likelihood:

[TABLE]

Furthermore, Eq. 8 can be further simplified under homoscedastic conditions:

[TABLE]

Hence, the probability density based loss function coincides with the well-known mean squared error (MSE).

2.3 Stochastic optimization for updating model parameters

Gradient based optimization is one of the most popular algorithms to optimize neural networks:

[TABLE]

where $\eta$ is generally known as the leraning rate that follows the Robbins-Monro conditions. The objective function stated in Eq. 11 indicates $\mathcal{O}\left(N\right)$ operations are required to compute $\mathcal{L}\left(\cdot\right)$ and $\nabla\mathcal{L}\left(\cdot\right)$ respectively, which may be computationally demanding for a large dataset. Therefore, stochastic gradient descent (SGD) is considered [14, 15, 32], where a random vector $g\left(\boldsymbol{x}\right)$ is defined to calculate the gradients [34]. With a restricted magnitude of stochastic gradients $\mathbb{E}||g\left(\boldsymbol{x}\right)||^{2}\leqslant N^{2}$ and a bounded variance $\mathbb{E}||g\left(\boldsymbol{x}\right)-\nabla\mathcal{L}\left(\boldsymbol{x}\right)||^{2}\leqslant\sigma^{2}$ [34], SGD can efficiently update model parameters by constructing a noisy natural gradient:

[TABLE]

In particular, adaptive moment estimation (ADAM) [35] that computes adaptive learning rates for $\boldsymbol{\omega}$ is adopted, and the corresponding $g\left(\boldsymbol{x}\right)$ takes the expression of:

[TABLE]

where $M_{t}$ and $V_{t}$ are estimates of the mean and variance of the gradients respectively. In ADAM, they are updated as follows [35]:

[TABLE]

It should be noted that the expresssion of Eq. 14 is an unbiased estimation of the exact gradient, and its calculation only depends on one data point.

3 Probabilistic modeling: a Bayesian approach

In this section, the aforementioned deterministic DL surrogate is enhanced to account for model uncertainties by the integration of Bayesian inference. Overall, Bayesian learning includes three steps: (1) establish prior beliefs about uncertain parameters; (2) compute the posterior distribution via Bayes’ rule; and (3) use the predictive distribution to determine a yet unobserved data point.

3.1 Prior representation: Bayesian hierarchical modelling

To begin with, let $\mathcal{U}_{\mathcal{E}}$ and $\mathcal{U}_{\mathcal{A}}$ represent the epistemic uncertainty and aleatory uncertainty respectively [8]. Prior information of $\mathcal{U}_{\mathcal{E}}$ and $\mathcal{U}_{\mathcal{A}}$ is initially encapsulated in a probability distribution function form. For the epistemic uncertainty, prior distributions are imposed on model parameters $\boldsymbol{\omega}$ [16, 17, 18]:

[TABLE]

Practical applications imply that the prior distribution $p\left(\boldsymbol{\omega}\right)$ should not be too restrictive on account of the limited prior information about $\boldsymbol{\omega}$ [20]. For this reason, hierarchical Bayesian modeling (HBM) method, which introduces a vector of hyperparameters $\boldsymbol{\eta}=[\eta_{1},\eta_{2},\dots,\eta_{n}]$ to the prior distribution, is employed to reduce subjective information induced undue influence on $p\left(\boldsymbol{\omega}\right)$ [9]. Consequently, the marginal prior can be obtained by integrating out $\boldsymbol{\eta}$ through the sum rule:

[TABLE]

where the joint probability distribution can be further expressed as a product of a set of conditional distributions via applying the product rule:

[TABLE]

For probabilistic modeling, model parameters in each layer of a BDL model are often assumed to follow a factorized multivariate Gaussian distribution:

[TABLE]

At the first hierarchy stage, let $\boldsymbol{\mu}_{\boldsymbol{\omega}}=\boldsymbol{0}$ and $\tau_{\boldsymbol{\omega}}$ be a Gamma random variable for instance. Hence, HBM breaks the prior distribution down to:

[TABLE]

with $\alpha_{\tau}$ and $\beta_{\tau}$ denoting the shape parameter and rate parameter of the Gamma distribution respectively. Using Eq. 17 and Eq. 18, the prior distribution can be reformulated as:

[TABLE]

In Eq. 21, $\mathcal{S}t\left(\cdot\right)$ characterizes the student’s t-distribution which is capable of providing heavier tails than Gaussian distribution.

Remark* (1).*

HBM grants a more impartial prior distribution by allowing the data to speak for itself [9], and it admits a more general modeling framework where the hierarchical prior becomes direct prior when the hyperparameters are modeled by a Dirac delta function (e.g. using $\delta\left(x-\tau_{\boldsymbol{\omega}}\right)$ to describe the precision term in Eq. 19). In addition, HBM offers the flexibility to work with a wide range of probability distributions, and even directly provides an analytical solution for some of the most popular choices such as Laplace, Gaussian, and student’s t-distribution [20].

On the other hand, homoscedastic noise $\boldsymbol{\epsilon}$ that is independent of the input data $\mathbb{V}ar[\boldsymbol{\epsilon}|\left(\boldsymbol{x}_{i},\boldsymbol{y}_{i}\right)]=\boldsymbol{\sigma}_{\boldsymbol{\epsilon}}^{2}\,\forall\,\left(\boldsymbol{x}_{i},\boldsymbol{y}_{i}\right)\in\boldsymbol{\mathcal{D}}$ is added to the output in consideration of the aleatory uncertainty $\mathcal{U}_{\mathcal{A}}$ which cannot be explained away by accepting more samples $\{\boldsymbol{x},\boldsymbol{y}\}$. Besides, additive noise term $\boldsymbol{\epsilon}$ guarantees a tractable likelihood for the probabilistic model $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ where $\boldsymbol{\epsilon}$ is most commonly modeled as a Gaussian process:

[TABLE]

Let $\boldsymbol{\mu}_{\boldsymbol{\epsilon}}=\boldsymbol{0}$ and $\boldsymbol{\sigma}_{\boldsymbol{\epsilon}}$ be a constant, the output vector hereby follows:

[TABLE]

A graphical model representation for the aforestated hierarchical prior as well as an illustration of the BDL model is given in Fig. 2.

3.2 Posterior approximation: variational inference

After defining a hierarchical prior distribution for the proposed probabilistic model $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$, the next step is to infer the posterior distribution, which reflects the updated parameter information. In the Bayesian formalism, the joint posterior distribution $p\left(\boldsymbol{\omega},\boldsymbol{\eta}\,|\,\boldsymbol{\mathcal{D}}\right)$ is calculated by:

[TABLE]

The marginal posterior distribution $p\left(\boldsymbol{\omega}\,|\,\boldsymbol{\mathcal{D}}\right)$ can be further determined by integrating out the joint posterior distribution, and the denominator of Eq. 24 is often referred to as the model evidence that takes the form of:

[TABLE]

In most cases, estimation of Eq. 25 it is computationally intractable as numerical integration requires a considerable number of samples if the parameter space $\mathcal{W}$ is very high [18]. To overcome this integration problem, variational inference (VI) is adopted so that Bayesian inference can proceed efficiently [26, 28].

Remark* (2).*

Different from the method of maximum a posteriori (MAP) which captures the mode of a posterior distribution [18], the objective for the posterior approximation at this place is to find a computationally efficient replacement of the true posterior distribution, so that numerical samples can be easily accessible in the predictive analysis.

3.2.1 Objective function: evidence lower bound

In VI, a family of proxy distributions parameterized by $\boldsymbol{\xi}$ is posited to approximate the true posterior distribution:

[TABLE]

VI attempts to make $q_{\boldsymbol{\xi}}\left(\boldsymbol{\omega}\right)$ looks as close as possible to $p\left(\boldsymbol{\omega}|\boldsymbol{\mathcal{D}}\right)$ via refining $\boldsymbol{\xi}$, and one typical interpretation of the closeness between two probability distributions is the Kullback-Leibler (KL) divergence [27, 28]. Therefore, VI casts the approximation problem in an optimization form, where the objective function can be expressed as:

[TABLE]

Instead of minimizing the KL divergence, we can equivalently maximize the evidence lower bound (ELBO) $\mathcal{L}\left(\boldsymbol{\omega}\right)$ [28]:

[TABLE]

The first conditional log-likelihood term in Eq. 28 is usually referred as to the data term [28, 29]. It compels the posterior distribution to explain data $\boldsymbol{\mathcal{D}}$ by maximizing the expected log-likelihood. Mini-batch optimization method is implemented to efficiently offer an unbiased stochastic estimator of the log-likelihood:

[TABLE]

where $M$ is a subset of $N$. Noticeably, besides accelerating the computational process, mini-batch optimization owns a higher model updating frequency that allows for a more robust convergence, and hence increases the chance of avoiding local minimum [34, 35]. Meanwhile, mean field variational inference (MFVI) method [26, 28, 29] is adopted to control the computational complexity of the second term in Eq. 28:

[TABLE]

The variational distribution is represented by a layer-wise factorized distribution where each factor is determined by its own variational parameter:

[TABLE]

Substituting Eq. 29 and Eq. 30 back to Eq. 28, the objective function of ELBO can be rewritten into:

[TABLE]

where the iteration of variational distribution for model parameters terminates when the convergence criteria is satisfied.

3.2.2 Gradients computation: stochastic gradient variational Bayes

Among the many techniques developed for solving optimization problems, gradient-based optimization method reliably tackles the EBLO maximization problem stated in Eq. 28 in an efficient manner [31]. For the sake of brevity, let:

[TABLE]

Using the log-derivative trick [28], the objective function $\mathcal{L}\left(\boldsymbol{\omega},\boldsymbol{\xi}\right)$ can be differentiated with respect to variational parameters $\boldsymbol{\xi}$:

[TABLE]

To quickly estimate numerical integrations, we can write Eq. 34 in its expectatio form and use Monte Carlo method to compute the stochastic gradients:

[TABLE]

However, it is observed that crude MC estimator for $\nabla_{\boldsymbol{\xi}}\mathcal{L}\left(\boldsymbol{\omega},\boldsymbol{\xi}\right)$ usually induces large variance [30, 31]. For this reason, stochastic gradient variational Bayes (SGVB) method is embraced to reduce the estimations’ variance [30]. Simply put, SGVB introduces an auxiliary variable $\boldsymbol{\epsilon}$ to the proxy distribution:

[TABLE]

where conditional probability density function $q_{\boldsymbol{\xi}}\left(\boldsymbol{\omega}|\boldsymbol{\epsilon}\right)$ is formally defined as a Dirac delta function:

[TABLE]

and $g\left(\boldsymbol{\xi},\boldsymbol{\epsilon}\right)$ is a differentiable transformation function that connects $\boldsymbol{\omega}$ and $\boldsymbol{\epsilon}$:

[TABLE]

For instance, a simple choice for $p\left(\boldsymbol{\epsilon}\right)$ is isotropic Gaussian distribution $\boldsymbol{\epsilon}\overset{i.i.d.}{\sim}p\left(\boldsymbol{\epsilon}\right)=\mathcal{N}\,\left(\boldsymbol{0},\,\boldsymbol{I}\right)$, and the reparameterization can be achieved though $\boldsymbol{\omega}=\boldsymbol{\mu}_{\boldsymbol{\omega}}+\boldsymbol{\sigma}_{\boldsymbol{\omega}}\odot\boldsymbol{\epsilon}$. Therefore, substituting Eq. 36 and Eq. 37 into Eq. 35, the pathwise estimator can be expressed as [30]:

[TABLE]

Combining Eq. 33 and Eq. 39, the final Monte Carlo estimator for the gradients can be written as:

[TABLE]

Now, the variance of stochastic gradients can be effectively reduced by magnitude of orders using this reparameterized estimator [30], and the VI-based optimization problem can be efficiently solved by the stochastic gradient descent algorithm mentioned in the previous section.

3.3 Predictive evaluation: Monte Carlo sampling

The last but the most important step of Bayesian computation concerns making predictions for new data samples $(\boldsymbol{x^{*}},\boldsymbol{y^{*}})$, where the predictive distribution can be expressed as:

[TABLE]

The optimized proxy posterior $q_{\boldsymbol{\xi}}\left(\boldsymbol{\omega}\right)$, which is obtained by solving the ELBO optimization problem, will take the place of the true posterior distribution $p\left(\boldsymbol{\omega}|\boldsymbol{\mathcal{D}}\right)$:

[TABLE]

In the same vein, the predictive integral is numerically achieved by drawing random samples from the proxy distribution. An unbiased estimator is given:

[TABLE]

For the purpose of uncertainty representation, it is of great importance to compute statistical moments of $\boldsymbol{y^{*}}$, such as mean:

[TABLE]

and variance:

[TABLE]

Because $\hat{\boldsymbol{y}}^{*}_{mean}$ and $\hat{\boldsymbol{y}}^{*}_{var}$ are essential elements for constructing the acquisition function which balances the exploration and exploitation in the context of Bayesian optimization [10, 13].

4 Numerical examples and results

4.1 Example 1: nonlinear regression

The first example considers a nonlinear function, which is commonly used as a testing problem to assess the accuracy of a regression model [13, 14, 36]. Mathematically, it is written as:

[TABLE]

To identify common features and differences between proposed model and current approaches, the regression problem is numerically solved using four different surrogate modeling methods: polynomial regression (PR); support vector machine (SVM); neural networks (NN); and Bayesian deep learning (BDL). First, a polynomial $\hat{f}\left(x\right)=\beta_{0}+\beta_{1}x+\dots+\beta_{n}x^{n}$ of degree $n=11$ is defined to fit the symmetric function [32]. The method of least squares is applied to find the best linear unbiased estimator (BLUE) of the regression coefficient vector $\boldsymbol{\beta}$ by minimizing the sum of squared errors. Secondly, an SVM regression model with a Gaussian kernel function $G\left(x_{i},x_{j}\right)=\exp\left(-\gamma||x_{i}-x_{j}||^{2}\right)$ is implemented to build a mapping between $x$ and $y$ [9, 10]. The default value for the kernel coefficient $\gamma$ is $1$, and sequential minimal optimization (SMO) algorithm is utilized to update the coefficients where Karush-Kuhn-Tucker (KKT) violation $\epsilon=0.0001$ is specified as the convergence criterion [37]. Thirdly, a feedforward neural network with one hidden layer that has $20$ neurons is built to learn the nonlinear transformation [14, 32]. Hyperbolic tangent function is adopted as the activation function for the hidden layer since its derivatives are steeper than sigmoid function. For the output layer, a straight line function that outputs the weighted sum from hidden neurons is used. Stochastic gradient descent is performed for parameter optimization [34], where the learning rate is fixed as a constant $\eta=0.0001$ and the default epoch setting is $100000$. Lastly, a Bayesian surrogate that has the same network configuration is examined. To account for the model uncertainty, a normal prior $\mathcal{N}\left(0,0.1\right)$ is directly imposed on the model parameters.

To train these models, we use the pseudorandom number generator to simulate a training dataset consisting of $30$ samples that are uniformly distributed in the interval $\left(-10,10\right)$. A Gaussian noise determined by $\epsilon_{i}\sim\mathcal{N}\left(0,\sigma\right)$ is added to each sample to make the problem more realistic [8]. Fig. 3 visualizes the fitted regression model via different approaches. It should be noted that the mean value of the predictive distribution is selected as the model estimation in the case of BDL. Obviously, BDL improves the generalization performance and mitigates the overfitting issue, which is encountered in NN modeling. Meanwhile, BDL is capable of characterizing the model uncertainty associated with the prediction in addition to achieving an equivalently accurate regression result compared to other methods. Table 1 and Table 2 summarize the coefficient of determination ($R^{2}$) and the root mean squared error (RMSE) for different surrogates. According to the results, BDL is more resistant to noisy data since increasing the random noise level deteriorates the effectiveness and quality of other three surrogates in a much more clear way.

4.2 Example 2: binary classification

To evaluate the classification performance of our proposed surrogate model, the second example applies the BDL to a synthetic dataset that holds a two-dimensional swirl pattern. As shown in Fig. 4, the synthetic dataset exhibits two intuitively separable manifolds, where each manifold resembles a crescent moon [9].

A BDL surrogate that is arranged in a $2\times 5\times 5\times 2$ form is developed as the neural network classifier. Specifically, two hidden layers are configured with the hyperbolic tangent activation function and softmax function $\sigma\left(\boldsymbol{x}_{i}\right)=\frac{e^{\boldsymbol{x}_{i}}}{\sum_{j=1}^{J}e^{\boldsymbol{x}_{j}}}$ is implemented to represent the categorical distribution for the outputs by computing a probability row vector where the sum of the row is $1$ [14, 15]. To understand the effects of different priors on the classification performance, we have considered three direct priors: Laplace $\mathcal{L}\left(0,1\right)$, Gaussian $\mathcal{N}\left(0,1\right)$, and Cauchy $\mathcal{C}\left(1,1\right)$. We further come up with three more hyper priors by fixing the location parameter $\alpha$ of the aforementioned probability distributions along with treating their scale parameter $\beta$ as a random variable, which can be described using an Inverse-Gamma distribution $\mathcal{IG}\left(1,1\right)$. The basic probability distribution functions are given as:

[TABLE]

Additionally, we conduct two trials to study the effects of noise on our neural network classifier. In the first trial, a BDL model is developed using $900$ samples, where each sample is contaminated by a Gaussian noise generated from $\epsilon_{i}\sim\mathcal{N}\left(0,0.1\right)$. In the second trial, $1200$ samples are utilized to build $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ as the external noise is amplified to $\epsilon_{i}\sim\mathcal{N}\left(0,0.3\right)$. Following the $70/30$ rule [14, 15, 32], the whole dataset $\mathcal{D}$ is divided into the training set $\mathcal{D}_{t}$ and the validation set $\mathcal{D}_{v}$, respectively. A first-order gradient-based optimization method, ADAM [35], is adopted to update model parameters, where the learning rate $\eta=0.001$, the exponential decay rates for the first/second moment estimates $\beta_{1}$ and $\beta_{2}$ are $0.9$ and $0.999$, respectively. It should be addressed that the reparametrization trick mentioned in Section 3.2 is automatically embedded by means of taking the derivatives of the objective function with respect to the variational parameters [30]. Here, the proxy distribution takes a Gaussian form, which indicates the variational posterior distribution is parameterized with two parameters, mean and standard deviation. To accelerate the training process, Mini-batch optimization method is used [14, 28], and the batch size is set to $30$. The stop criteria epoch number is $50000$.

Fig. 5 provides a graphic illustration of the classification results. According to these results, BDL model becomes less confident about its predictions when validation samples are more near the true separation trajectory. It is because even small noise can distort the original manifold in a severe way [6]. However, the proposed hyper priors are able to provide better predictions especially in the second trial where the addictive noise is stronger. This is credited to the nature mechanism of Bayesian hierarchical modeling, which relaxes the prior constraints by encoding prior belief using a series of hyperparameter values instead of fixed constants [9]. Lastly, Fig. 6 reveals the variational posterior distribution of weights and bias in the first hidden layer. For the previous proposed priors, zero centered Laplace prior is equivalent to the $L1$ regularization and Gaussian prior is identical to the $L2$ regularization [20, 25]. In Fig. 6, results of $\mathcal{L}\left(0,1\right)$ is approximately sparse signal and $p\left(\boldsymbol{\omega}\right)$ of $\mathcal{N}\left(0,1\right)$ is not centering around zero, which aligns with the properties of $L1$ and $L2$ regularization respectively [32].

4.3 Example 3: structural analysis of a geometrically nonlinear membrane

This example addresses the computational cost issue of using finite element (FE) model in the structural analysis with uncertain inputs [38]. The target structure is a geometrically nonlinear membrane that is clamped at four edges [39, 40]. Fig. 7. (a.1) gives a sketch of the objective domain $\Omega=[0,l]\times[0,b]\subset\mathbb{R}^{2}$, where uniformly distributed pressure loads are applied on the upper surface. We are interested in using the BDL based surrogate $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ to approximate the nonlinear mechanism between uncertain structure parameters $\boldsymbol{x}$ and random responses $\boldsymbol{y}$.

4.3.1 Membrane example: nonlinear analysis and surrogate modeling

Uncertainty analysis. Vector $\boldsymbol{x}$ covers geometric uncertainties, where $x_{1}=l$, $x_{2}=b$, and $x_{3}=t$ are the length, breadth, and thickness of target membrane, as well as material uncertainties, with $x_{4}=E$ and $x_{5}=\upsilon$ denoting the elastic modulus and Poisson’s ratio, respectively. Table 3 gives a systematic summary of statistical properties of $\boldsymbol{x}$. The quantities of interest $\boldsymbol{y}$ are z-direction displacements $w$ at locations of $p_{1}-\left(1,0.5\right)$, $p_{2}-\left(0.6,0.2\right)$, and $p_{3}-\left(1.6,0.8\right)$. The load-displacement relationship is no longer a deterministic curve due to the input randomness (See Fig. 7. (a.2)).

Nonlinear FE model. Because of the geometric nonlinearity, the inplain strain is partitioned into two parts $\boldsymbol{\epsilon}=\boldsymbol{\epsilon}^{l}+\boldsymbol{\epsilon}^{non}$, where $\boldsymbol{\epsilon}^{l}$ describes the linear strain and $\boldsymbol{\epsilon}^{non}$ represents the nonlinear strain term:

[TABLE]

The solution $\boldsymbol{d}=[u,v,w]$ of these nonlinear equilibrium equations are obtained by the Newton-Raphson (NR) method [39]. The iterative process terminates when the unbalanced force residual is smaller than the tolerance $\epsilon=0.0001$ or the NR algorithm reaches the default maximum iteration $n=100$. The force and displacement vector is initialized to zero and the increment loads $\Delta_{p}=p/n$ where $p=100$ and $n=400$. It should be noted that the thickness of the membrane is comparatively small in relation to other two dimensions. It is therefore the FE model $\mathcal{M}\left(\cdot\right)$ can be built by 200 $(20\,\text{in x axis}\times 10\,\text{in y axis})$ four node (Q4) quadrilateral elements and large deformation theory is adopted [39].

Surrogate model. We use the proposed BDL approach to provide a $\mathbb{R}^{5}\rightarrow\mathbb{R}^{3}$ transformation. The network architecture has three hidden layers $30\times 15\times 10$ besides the input and output layer [33] and probability distributions that account for model uncertainties are specified in a layer by layer fashion. To investigate the efficacy of different prior, a direct zero-mean Gaussian $\mathcal{N}\left(0,1\right)$ and a hierarchical prior $\mathcal{N}\left(0,\mathcal{IG}\left(1,1\right)\right)$ have been tested. Furthermore, training datasets of size $64$, $128$, $256$, and $512$ have been considered for the purpose of identifying the influence from the amount of data on the accuracy of model predictions. For the posterior approximation, ADAM is adopted [35], where the initial learning rate $\eta$ is $0.005$. Notably, $\eta$ decays every $100$ epochs by multiplying a constant rate of $0.75$ and the epoch number is $1000$. The batch size for the subsampling procedure of all trials is set to $16$. In the variational inference stage, $200$ numerical samples are employed to estimate the lower bound.

4.3.2 Results

In Fig. 7, the optimization results imply the predictive distribution computed by Eq. 43 shrinks rapidly as the number of training sample increases. $128$ samples can give a narrow-band distribution of $w\left(p_{1}\right)$, indicating the trained BDL model becomes sufficiently reliable as most model uncertainties have been explained away by data. Moreover, we compared the evolution process of the predictive distribution $p\left(w_{p_{1}}\right)$ via different priors. In both trials, $128$ training samples are used, and the same validation sample is randomly chosen where $\boldsymbol{x}=[2.0117,1.0157,0.0019,213.1180,0.3018]$ and $\boldsymbol{y}=[-2.3428,-1.1501,-1.0536]$. It is found that the predictive distribution via hyper prior takes more epochs to shrink (See Fig. 7. (c.1) and (c.2)). The intuitive explanation is $\mathcal{N}\left(0,\mathcal{IG}\left(1,1\right)\right)$ has a larger initial parameter space than $\mathcal{N}\left(0,1\right)$. Despite the difference, both BDL surrogates provide a reliable input-output mapping, and the coefficient of determination is summarized in Table 4. To check the generalizability of proposed surrogates, $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ is further applied to uncertainty analysis. First, $1000$ samples were used to train the network. Next, another $1\times 10^{5}$ samples are exploited to develop the distribution of displacements at $p_{1}$, $p_{2}$, and $p_{3}$. Fig. 8 presents the UQ results, where BDL based surrogates accurately propagate uncertainties to the response distributions.

4.4 Example 4: prediction of wind pressure

Obtaining detailed data of wind-induced pressure coefficients on building surfaces is of great practical importance in the design of high-rise buildings. However, wind tunnel test results are limited and may be contaminated through different sources. In this example, the proposed BDL model is applied to predict the mean and root-mean-square (RMS) pressure coefficients using limited experiment data.

4.4.1 Wind pressure database and predictive model

Wind tunnel data. The aerodynamic database considered for this example is developed by the Tokyo Polytechnic University (TPU) [41]. For the wind tunnel experiment, a $1:400$ scale rigid model was built to represent the target tall building of dimension $200m\times 40m\times 40m$ and a power law exponent of $1/4$ was used for the description of the mean wind speed. A total of $500$ pressure taps were used to collect data at a sampling frequency of $1000\,Hz$ for a sample period of $32.768\,s$. Hourly average wind speed was $11.1438\,m/s$ and wind attacking angle was $0^{\circ}$, indicating wind direction is perpendicular to the front face of the experiment model. The wind pressure of our interest is characterized by a dimensionless number $\boldsymbol{C_{p}}\left(\boldsymbol{x}\right)=\frac{p_{\boldsymbol{x}}-p_{\infty}}{\rho U_{0}^{2}/2}$ that is known as the pressure coefficient. $p_{\infty}$ is the static pressure at freestream, $\rho$ is the air density and $U_{0}$ is the mean wind speed at the reference height. The predictive quantities are:

[TABLE]

To demonstrate the efficacy of the BDL surrogate in dealing with small datasets, Biharmonic spline interpolation is performed based on the measured pressure data $N_{old}=500$. Specifically, the are $250$ width interpolation points and $750$ height interpolation points in each building face. As a result, the surface pressure fields are described by a total of $N_{new}=250\times 750\times 4=750000$ synchronous pressure points. However, we only use as much as $1\%$ of the total data to train the Bayesian model.

Prediction model. The Cartesian coordinates $\boldsymbol{x}\in\mathbb{R}^{2}$ are selected input variables of the BDL model, and the output is a scalar either $y=C_{p}^{mean}\left(\boldsymbol{x}\right)$ or $y=C_{p}^{rms}\left(\boldsymbol{x}\right)$. After extensive hyperparameters and network architectures search, it was found that BDL with network configuration of $2\times 15\times 10\times 1$ and $2\times 30\times 15\times 1$ provide superior performance in predicting $C_{p}^{mean}$ and $C_{p}^{rms}$, respectively. Hyperbolic tangent function $tanh\left(x\right)=\frac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$ is adopted as the activation function since it produces steep derivatives, and ADMA optimizer [35] is implemented with a learning rate initialized to $0.03$, which follows a step decay to prevent optimizing parameters chaotically. The annealing strategy is adopted to improve the stochastic gradients where the anneal rate is fixed to $0.75$. $1000$ samples are applied to get a reasonable estimation for the test log likelihood during the training phase. Epoch number is set as $1000$ and the testing frequency is $10$. To validate the effectiveness of HBM based prior, a direct Gaussian prior $\mathcal{N}\left(0,1\right)$ and a hyper prior $\mathcal{N}\left(0,\mathcal{IG}\left(1,1\right)\right)$ are examined.

4.4.2 Results

Fig. 9 summarizes the $\boldsymbol{C_{p}}^{mean}$ predictions and Fig. 10 provides the predictive results of $\boldsymbol{C_{p}}^{rms}$. The ground truth is numerically obtained via the interpolation and extrapolation of the wind tunnel data $\boldsymbol{C_{p}}$; the predictions are defined as the mean value of the predictive distribution $\hat{\boldsymbol{C_{p}}}$; the relative error measures the predictive difference between $\boldsymbol{C_{p}}$ and $\hat{\boldsymbol{C_{p}}}$; and model uncertainties are evaluated using the variance of $\hat{\boldsymbol{C_{p}}}$. The results reveal that hyperprior $\mathcal{N}\left(0,\mathcal{IG}\left(1,1\right)\right)$ not only allows $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ to be better fitted but also makes $\mathcal{M}^{{}^{\prime}}\left(\cdot\right)$ less sensitive to the complex nature of the noise embedded in the experimental data. It can be seen from the limits of the colorbar that the magnitude of relative errors is reduced for those hyperprior trials. From a percent error perspective, most predictions are lying within the $97\%$ confidence interval and the maximum errors are all less than $10\%$, most of which are mainly gathered around the boundary due to the extrapolation algorithm that is performed at the data preprocessing stage. In Fig. 9, uncertainties results confirm the propagation of uncertainty resulting from the extrapolation process and successfully identify the areas that are prone to. In general, it seems hyperprior improves the model performance more in the case of predicting $\boldsymbol{C_{p}}^{rms}$ than $\boldsymbol{C_{p}}^{mean}$, which is reasonable as the second sample moment magnifies the differences between predicted values and observed values, causing the need of more degrees of freedom to explain the variation. Table 5 summarizes the model performance using RMSE. To illustrate the robustness of the variational inference method in terms of approximating intractable posterior distributions, Fig. 11 shows the optimization process of computing the variational posterior. It can be observed that VI is capable of capturing the main characteristics within the first $30$ epochs, and it is applicable to various situations where the true posterior distribution takes different kinds of forms.

5 Concluding Remarks

This paper presents a probabilistic modeling approach for learning hidden relationships from limited and noisy data using Bayesian deep learning (BDL) with hierarchical prior. The proposed surrogate rigorously accounts for the model uncertainties by means of imposing prior distributions on model parameters. Meanwhile, it effectively propagates the preassigned prior belief to the prediction quantities. In summary, the following conclusions are drawn:

• (1)

Bayesian inference has been successfully integrated into the current deterministic deep learning framework. Consequently, the proposed model is able to analyze uncertainties associated with model predictions and help stakeholders make a more informed decision by providing a confidence level for the predictive estimation.

• (2)

The hypothesis of using hierarchical Bayesian modeling to describe prior distributions of model parameters is tested. In both classification and regression problems, superior performances can be achieved utilizing hyper priors, especially when the training data is seriously contaminated. Moreover, probabilistic surrogate with hyper prior tends to have an improved learning ability from a small dataset.

• (3)

Intractable posterior distributions that risen from multidimensional integrals step of Bayesian analysis has been addressed by the state-of-the-art variational inference method. Compared to some advanced sampling-based methods, variational inference method offers a higher scalability by tackling the model learning problem in an objective-equivalent-transformed and gradients-effective-computed optimization form.

• (4)

The examples provided have demonstrated the applicability of the proposed modeling scheme to both classification and regression tasks involving complex systems. Especially in the membrane example, BDL is capable of providing an accurate description of the highly nonlinear mapping between different design variables and various structural performance indicators and produces virtually identical uncertainty quantification results as conventional Monte Carlo method. Furthermore in the wind field prediction example, BDL model is trained and tested using very limited wind tunnel data, and it is shown that the probabilistic model is not only able to effectively recover the entire mapping of the mean and root-mean-square pressure fields with high precision using as small as $1\%$ of the data but also quantifies the uncertainty level at every single point in the prediction domain, serving as a reliable surrogate for learning complex field distribution.

To improve the model performance, it is envisaged that the combination of the BDL model with information from underlying physics can not only further accelerate the training of neural networks but also holds the promise of interpreting the learning process.

6 Acknowledgments

This work was supported by the National Science Foundation (NSF) under Grant No. 1520817 and No. 1612843. This support is gratefully acknowledged.

7 Reference