Uncertainty quantification and reliability analysis by an adaptive sparse Bayesian inference based PCE model

Abstract

An adaptive Bayesian polynomial chaos expansion (BPCE) is developed in this paper for uncertainty quantification (UQ) and reliability analysis. The sparsity in the PCE model is developed using automatic relevance determination (ARD) and the PCE coefficients are computed using the variational Bayesian (VB) inference. Further, Sobol sequence is utilized to evaluate a response quantity sequentially. Finally, leave one out (LOO) error is used to obtain the adaptive BPCE model. UQ and reliability analysis are performed of some numerical examples by the adaptive BPCE model. It is found that the optimal number of model evaluations and the optimal PCE degree are suitably selected simultaneously for a problem by the adaptive BPCE model. A highly accurate result is predicted by the proposed approach using very few model evaluation. Further, highly sparse PCE models are obtained by the ARD approach for most of the numerical examples. Additionally, distribution parameters of the predicted response quantity are also obtained by the VB inference, which are used to compute the confidence interval of the predicted response quantities.

Appendices

Appendix A: Formulation of VLB by mean field theory

Recall the VLB as given in Eq. (10):

$$\begin{aligned}&{\mathcal {L}}\left[ q\left( \varTheta \right) \right] \nonumber \\&\quad = \int q\left( \varTheta \right) \log \frac{p\left( \varTheta ,Y\right) }{q\left( \varTheta \right) } {\text {d}}\varTheta \end{aligned}$$

(51)

$$\begin{aligned}&\quad = \int q\left( \varTheta \right) \left[ \log p\left( \varTheta ,Y\right) -\log q\left( \varTheta \right) \right] {\text {d}}\varTheta \end{aligned}$$

(52)

$$\begin{aligned}&\quad = \int \prod _{i=1}^{n_v} q\left( \varTheta _i\right) \left[ \log p\left( \varTheta ,Y\right) -\sum _{i=1}^{n_v} \log q\left( \varTheta _i\right) \right] {\text {d}}\varTheta \end{aligned}$$

(53)

$$\begin{aligned}&\quad = \int q\left( \varTheta _i\right) \left[ \int \log p\left( \varTheta ,Y\right) \prod _{j\ne i}q\left( \varTheta _j\right) \right] {\text {d}}\varTheta _i \nonumber \\&\qquad - \int q\left( \varTheta _i\right) \log q\left( \varTheta _i\right) {\text {d}}\varTheta _i + {{\,\mathrm{const}\,}}\end{aligned}$$

(54)

$$\begin{aligned}&\quad = \int q\left( \varTheta _i\right) \log {\tilde{p}} \left( \varTheta _i,Y\right) {\text {d}}\varTheta _i \nonumber \\&\qquad - \int q\left( \varTheta _i\right) \log q\left( \varTheta _i\right) {\text {d}}\varTheta _i + {{\,\mathrm{const}\,}}\end{aligned}$$

(55)

$$\begin{aligned}&\quad = -{{\,\mathrm{KL}\,}}\left( q\left( \varTheta _i\right) \parallel {\tilde{p}}\left( \varTheta _i,Y\right) \right) + {{\,\mathrm{const}\,}}. \end{aligned}$$

(56)

A new probability distribution ${\tilde{p}}\left( \varTheta _i,Y\right)$ is defined in Eq. (55): it is also represented as follows [9]:

$$\begin{aligned} \log {\tilde{p}}\left( \varTheta _i,Y\right) = {\mathbb {E}}_{j\ne i} \left[ \log p \left( \varTheta ,Y\right) \right] , \end{aligned}$$

(57)

where ${\mathbb {E}}\left[ \bullet \right]$ represents the expectation with respect to the j-th parameters such that $j\ne i$. It is evident from Eq. (51) that the VLB is maximum when the KL divergence is minimum: the KL divergence is minimum when $q\left( \varTheta _i\right) = {\tilde{p}}\left( \varTheta _i,Y\right)$. Hence, the minimum KL divergence for the i-th parameter occurs at:

$$\begin{aligned} \log q\left( \varTheta _i\right) = {\mathbb {E}}_{j\ne i} \left[ \log p\left( \varTheta ,Y\right) \right] . \end{aligned}$$

(58)

Appendix B: Distributions for likelihood function and prior

The posterior is best estimated considering the mixed distribution [9]. Hence, the likelihood function is inferred by the normal distribution:

$$\begin{aligned} p\left( {\mathcal {Y}}|\varPsi ,a,\varsigma \right)= & {} \prod _{i=1}^{N} {\mathcal {N}}\left( {\mathcal {Y}}_i|\varPsi _i a,\varsigma ^{-1}\right) \end{aligned}$$

(59)

$$\begin{aligned}= & {} \left( \frac{\varsigma }{2\pi }\right) ^{\frac{N}{2}} \exp \left[ -\frac{\varsigma }{2}\sum _{i=1}^{N}\left( {\mathcal {Y}}_i-\varPsi _i a\right) ^2\right] , \end{aligned}$$

(60)

where ${\mathcal {N}}\left( \bullet \right)$ is referred as the normal distribution. To maintain the conjugacy [9], the prior is inferred by a conjugate normal-gamma distribution [25]:

$$\begin{aligned} p\left( a,\varsigma |\omega \right)= & {} {\mathcal {N}}\left( a|0,\left( \varsigma \omega \right) ^{-1}\right) {{\,\mathrm{Gam}\,}}\left( \varsigma |A_0,B_0\right) \end{aligned}$$

(61)

$$\begin{aligned}= & {} \left( 2\pi \right) ^{-\frac{n}{2}} \left| \varvec{\omega }\right| ^\frac{1}{2} \frac{B_0^{A_0}}{\varGamma \left( A_0\right) } \varsigma ^{\frac{n}{2}+A_0-1}\nonumber \\&\exp \left[ -\frac{\varsigma }{2}\left( a^T\varvec{\omega } a+2B_0\right) \right] , \end{aligned}$$

(62)

where ${{\,\mathrm{Gam}\,}}\left( \bullet \right)$ is the gamma distribution and $A_0,B_0$ are the gamma distribution parameters for $\varsigma$. $\omega = \{\omega _1,\omega _2,\ldots ,\omega _{n}\}^T = {{\,\mathrm{diag}\,}}\left( \varvec{\omega }\right)$ is the hyper-prior. $\varvec{\omega } \in {\mathbb {R}}^{n \times n}$ is the sparse matrix having only the diagonal terms. Further, hyper-prior is inferred by a gamma distribution:

$$\begin{aligned} p\left( \omega \right)= & {} \prod _{j=1}^{n} {{\,\mathrm{Gam}\,}}\left( \omega _j|C_0,D_0\right) \end{aligned}$$

(63)

$$\begin{aligned}= & {} \prod _{j=1}^{n} \frac{D_0^{C_0}}{\varGamma \left( C_0\right) } \omega _j^{C_0-1} \exp \left( -D_0\omega _j\right) , \end{aligned}$$

(64)

where $C_0,D_0$ are the parameters of the gamma distribution for the hyper-prior.

Appendix C: Predictive distribution

The predictive distribution at the new samples $\varXi _{{{\text {pred}}}}$ can be obtained having the available informations ${\mathcal {D}}\in \left\{ \varXi ,{\mathcal {Y}} \right\}$ [9]. Marginalizing over the parameters, the predictive distribution is:

$$\begin{aligned}&p\left( {\mathcal {Y}}_{{{\text {pred}}}}|{ {\hat{\varPsi }} }_{{{\text {pred}}}},{\mathcal {D}} \right) \nonumber \\&\quad =\iiint { p\left( {\mathcal {Y}}_{{{\text {pred}}}}|{ {\hat{\varPsi }} }_{{{\text {pred}}}},a,\varsigma \right) p\left( a,\varsigma ,\omega |{\mathcal {D}} \right) {\text {d}} a {\text {d}}\varsigma {\text {d}}\omega } \end{aligned}$$

(65)

$$\begin{aligned}&\quad \approx \iiint { p\left( {\mathcal {Y}}_{{{\text {pred}}}}|{ {\hat{\varPsi }} }_{{{\text {pred}}}},a,\varsigma \right) q\left( a,\varsigma \right) q\left( \omega \right) {\text {d}} a {\text {d}}\varsigma d\omega } \end{aligned}$$

(66)

$$\begin{aligned}&\quad \iint {\mathcal {N}}\left( {\mathcal {Y}}_{{{\text {pred}}}}|{ {\hat{\varPsi }} }_{{{\text {pred}}}}{ a }_{ k },{ \varsigma }^{ -1 } \right) {\mathcal {N}}\nonumber \\&\qquad \left( a|{ a }_{ k },{ \varsigma }^{ -1 }{ h }_{ k } \right) {{\,\mathrm{Gam}\,}}\left( \varsigma |{ A }_{ k },{ B }_{ k } \right) {\text {d}} a {\text {d}}\varsigma \end{aligned}$$

(67)

$$\begin{aligned}&\quad =\int {\mathcal {N}}\left( {\mathcal {Y}}_{{{\text {pred}}}}|{ {\hat{\varPsi }} }_{{{\text {pred}}}}{ a }_{ k },{ \varsigma }^{ -1 }\right. \nonumber \\&\qquad \left. \left( 1+{ {\hat{\varPsi }} }_{{{\text {pred}}}}{ h }_{ k }{ {\hat{\varPsi }} }_{{{\text {pred}}}}^{ T } \right) \right) {{\,\mathrm{Gam}\,}}\left( \varsigma |{ A }_{ k },{ B }_{ k } \right) {\text {d}}\varsigma \end{aligned}$$

(68)

$$\begin{aligned}&\quad ={{\,\mathrm{St}\,}}\left( {\mathcal {Y}}_{{{\text {pred}}}}|\mu ,\lambda ,\nu \right) . \end{aligned}$$

(69)

It is seen from Eq. (65) that the predictive distribution is Student’s t distribution (denoted by ${{\,\mathrm{St}\,}}$) with the parameters $\mu$, $\lambda$ and $\nu$. After constructing the adaptive BPCE model, the distribution parameters are given by:

$$\begin{aligned}&\mu ={ {\hat{\varPsi }} }_{{{\text {pred}}}} {\hat{a}} \end{aligned}$$

(70)

$$\begin{aligned}&\lambda =\frac{ { A }_{ k } }{ { B }_{ k } } { \left( 1+{ {\hat{\varPsi }} }_{{{\text {pred}}}}{ h }_{ k }{ {\hat{\varPsi }} }_{{{\text {pred}}}}^T \right) }^{ -1 } \end{aligned}$$

(71)

$$\begin{aligned}&\nu =2{ A }_{ k }. \end{aligned}$$

(72)

The standard deviation of the predictive distribution is given by:

$$\begin{aligned} \sigma = \sqrt{{ \left( 1+{ {\hat{\varPsi }} }_{{{\text {pred}}}}{ h }_{ k }{ {\hat{\varPsi }} }_{{{\text {pred}}}}^{ T } \right) } \frac{B_k}{A_k -1}}. \end{aligned}$$

(73)