Global sensitivity analysis using polynomial chaos expansion enhanced Gaussian process regression method

Abstract

Global sensitivity analysis (GSA) is a commonly used approach to explore the contribution of input variables to the model output and identify the most important variables. However, performing GSA typically requires a large number of model evaluations, which can result in a heavy computational burden, particularly when the model is computationally expensive. To address this issue, an efficient Sobol index estimator is proposed in this paper using polynomial chaos expansion (PCE) enhanced Gaussian process regression (GPR) method, namely PCEGPR. The orthogonal polynomial functions of PCE method are incorporated into GPR surrogate model to construct the kernel function. An estimation scheme based on fixed-point iteration and leave-one-out cross-validation error is presented to determine the optimal parameters of PCEGPR method. The analytical expressions of main and total sensitivity indices are also derived by considering the posterior predictor and covariance of PCEGPR surrogate model. The effectiveness of the proposed estimator is demonstrated by four numerical examples.

Appendix

In the general case, the random noise $\upsilon$ with zero mean and $\sigma^{2}$ variance, i.e., $\upsilon \sim \boldsymbol{\mathcal{N}}\left( {0,\sigma^{2} } \right)$ is considered in the training data. Then, the covariance matrix in Eq. (20) can be written as $\Lambda = \Phi^{T} W\Phi + \sigma^{2}I$, where $I$ denotes the identical matrix. The maximum likelihood method is used to compute $\mu$ and $W$. The likelihood function is written as

$$p\left( {\boldsymbol{\mathcal{Y}}|\mu ,W} \right) = \left( {2\pi } \right)^{ - N/2} \left| \Lambda \right|^{ - 1/2} \exp \left( { - \frac{{{\varvec{\beta}}^{T} \Lambda^{ - 1} {\varvec{\beta}}}}{2}} \right)$$

(66)

Taking the logarithm operator of likelihood function, the log-likelihood function is given as

$$\boldsymbol{\mathcal{L}}\left( {\mu ,W} \right) = \ln p\left( {\boldsymbol{\mathcal{Y}}|\mu ,W} \right) = - \frac{N}{2}\ln \left( {2\pi } \right) - \frac{1}{2}\ln \left| \Lambda \right| - \frac{1}{2}{\varvec{\beta}}^{T} \Lambda^{ - 1} {\varvec{\beta}}$$

(67)

Taking the derivative of $\boldsymbol{\mathcal{L}}\left( {\mu ,W} \right)$ with respect to $\mu$ and setting it to zero, we have

$$\frac{{\partial \boldsymbol{\mathcal{L}}\left( {\mu ,W} \right)}}{\partial \mu } = \frac{{\partial \beta^{T} \Lambda^{ - 1} \beta }}{\partial \mu } = 2F^{T} \Lambda^{ - 1} \left( {\boldsymbol{\mathcal{Y}} - F\mu } \right) = 0$$

(68)

Then, the maximum likelihood estimation of $\mu$ is given as

$$\hat{\mu } = \left( {F^{T} \Lambda^{ - 1} F} \right)^{ - 1} F^{T} \Lambda^{ - 1} \boldsymbol{\mathcal{Y}}$$

(69)

The second term of $\boldsymbol{\mathcal{L}}\left( {\mu ,W} \right)$ in Eq. (67) can be written as

$$\ln \left| \Lambda \right| = \ln \left( {\left| {\sigma^{2} I} \right|\left| C \right|\left| W \right|} \right) = \ln \left| C \right| + N\ln \sigma^{2} + \ln \left| W \right|$$

(70)

where $C = \sigma^{ - 2} \Phi \Phi^{T} + W^{ - 1}$. According to Woodbury inversion [54], $\Lambda^{ - 1}$ is reformulated as

$$\Lambda^{ - 1} = \left( {\Phi^{T} W\Phi + \sigma^{2}I } \right)^{ - 1} = \sigma^{ - 2}I - \sigma^{ - 2} \Phi^{T} C^{ - 1} \Phi \sigma^{ - 2}$$

(71)

and the vector of PCE coefficients PCE in Eq. (25) is also deduced as

$${\varvec{\lambda}} = W\Phi \left( {\Phi^{T} W\Phi + \sigma^{2}I} \right)^{ - 1} {\varvec{\beta}} = \sigma^{ - 2} C^{ - 1} \Phi {\varvec{\beta}}$$

(72)

Based on Eqs. (71) and (72), the third term of $\boldsymbol{\mathcal{L}}\left( {\mu ,W} \right)$ is derived as [49, 54]

$${\varvec{\beta}}^{T} \Lambda^{ - 1} {\varvec{\beta}} = \sigma^{ - 2} \left( {{\varvec{\beta}} - \Phi^{T} {\varvec{\lambda}}} \right)^{T} \left( {{\varvec{\beta}} - \Phi^{T} {\varvec{\lambda}}} \right) + {\varvec{\lambda}}^{T} W^{ - 1} {\varvec{\lambda}}$$

(73)

Taking the derivative of $\boldsymbol{\mathcal{L}}\left( {\mu ,W} \right)$ with respect to $w_{{{\varvec{\alpha}}_{i} }}$ and forcing it to zero, the following equation is obtained

$$\frac{{\partial \boldsymbol{\mathcal{L}}\left( {\mu ,W} \right)}}{{\partial w_{{{\varvec{\alpha}}_{i} }} }} = \frac{{C_{{{\varvec{\alpha}}_{i} {\varvec{\alpha}}_{i} }}^{ - 1} - w_{{{\varvec{\alpha}}_{i} }}^{2} + \lambda_{{{\varvec{\alpha}}_{i} }}^{2} }}{{w_{{{\varvec{\alpha}}_{i} }}^{3} }} = 0$$

(74)

where $C_{{{\varvec{\alpha}}_{i} {\varvec{\alpha}}_{i} }}^{ - 1}$ is the ith diagonal element of $C^{ - 1}$. Then the estimation of $w_{{{\varvec{\alpha}}_{i} }}$ is given as

$$\hat{w}_{{{\varvec{\alpha}}_{i} }} = \left( {C_{{{\varvec{\alpha}}_{i} {\varvec{\alpha}}_{i} }}^{ - 1} + \lambda_{{{\varvec{\alpha}}_{i} }}^{2} } \right)^{1/2}$$

(75)

In Sect. 4.1, the random noise is not taken into consideration in the GPR modeling, i.e., $\sigma^{2} = 0$. Therefore, Eq. (69) is transformed as

$$\hat{\mu } = \left( {F^{T} R^{ - 1} F} \right)^{ - 1} F^{T} R^{ - 1} \boldsymbol{\mathcal{Y}}$$

(76)

and $C^{ - 1} = \sigma^{2} \left( {\Phi \Phi^{T} + \sigma^{2} W^{ - 1} } \right)^{ - 1}$ is equal to 0. Then, $\hat{w}_{{{\varvec{\alpha}}_{i} }} = \lambda_{{{\varvec{\alpha}}_{i} }}$ and

$$\hat{\user2{w}} = {\varvec{\lambda}} = W\Phi \left( {\Phi^{T} W\Phi } \right)^{ - 1} {\varvec{\beta}} = W\Phi {\varvec{\tau}}$$

(77)

To compute the values of $\hat{\mu }$ and $\hat{\user2{w}}$, the fixed-point iteration is used to solve Eqs. (76) and (77) in Sect. 4.1.