A non-intrusive reduced-order modeling method using polynomial chaos expansion

TL;DR

This paper introduces a non-intrusive reduced-order modeling approach combining POD and PCE to efficiently perform uncertainty quantification with reduced computational costs and high accuracy.

Contribution

It presents a novel non-intrusive method integrating POD and PCE for stochastic modeling, reducing costs while maintaining accuracy in UQ analysis.

Findings

Significantly reduces computational costs and storage requirements.

Achieves similar accuracy to full-PCE in statistical predictions.

Effectively predicts outputs using only a few snapshots.

Abstract

We propose a non-intrusive reduced-order modeling method based on proper orthogonal decomposition (POD) and polynomial chaos expansion (PCE) for stochastic representations in uncertainty quantification (UQ) analysis. Firstly, POD provides an optimally ordered basis from a set of selected full-order snapshots. Truncating this optimal basis, we construct a reduced-order model with undetermined coefficients. Then, PCE is utilized to approximate the coefficients of the truncated basis. In the proposed method, we construct a PCE using a non-intrusive regression-based method. Combined with the model reduction ability of POD, the proposed method efficiently provides stochastic representations in UQ analysis. To investigate the performance of the proposed method, we provide three numerical examples, i.e., a highly nonlinear analytical function with three uncertain parameters, two-dimensional (2D) heat-driven cavity flow with a stochastic boundary temperature, and 2D heat diffusion with stochastic conductivity. The results demonstrate that the proposed method significantly reduces the computational costs and storage requirements that arise due to high-dimensional physical and random spaces. While demonstrating a similar accuracy with that of the classical full-PCE in predicting statistical quantities. Furthermore, the proposed method reasonably predict the outputs of the full order model using only a few snapshots.

Figures

Tables

Table 1. Table 1: Effects of ε 𝜀 \varepsilon on L 𝐿 L and MAREs in mean and variance predictions for N s = 400 subscript 𝑁 𝑠 400 N_{s}=400 for the stochastic Ackley function.

Tolerance $\epsilon$		Size of selected POD basis $L$		MARE in mean			MARE in variance
				${\epsilon }_{r,\mu }^{MC}$	${\epsilon }_{r,\mu }^{FPCE}$		${\epsilon }_{r,var}^{MC}$	${\epsilon }_{r,var}^{FPCE}$
${10}^{-4}$		7		2.4695E-04	2.4117E-04		1.8000E-02	1.6200E-02
${10}^{-5}$		13		6.0376E-05	5.4439E-05		5.4000E-03	2.3000E-03
${10}^{-6}$		20		2.6024E-05	1.7988E-05		4.5000E-03	6.7358E-04
${10}^{-7}$		29		1.9915E-05	1.0626E-05		4.3000E-03	2.5562E-04
${10}^{-8}$		40		1.6268E-05	3.2476E-06		4.3000E-03	1.1670E-04
${10}^{-9}$		52		1.6278E-05	2.3119E-06		4.3000E-03	7.7865E-05

Table 2. Table 2: Effects of N s subscript 𝑁 𝑠 N_{s} on L 𝐿 L and MAREs in mean and variance predictions for ε = 10 − 8 𝜀 superscript 10 8 \varepsilon=10^{-8} for the stochastic Ackley function.

No. of snapshots		Size of selected POD basis		MARE in mean		MARE in variance
$N_s$		$L$		${\epsilon }_{r,\mu }^{FPCE}$		${\epsilon }_{r,var}^{FPCE}$
50		31		6.2769E-04		3.1000E-02
100		37		2.0937E-04		4.0000E-03
200		39		4.0501E-05		9.0066E-03
400		40		3.2476E-06		1.1670E-04
600		41		2.6858E-06		9.7453E-05

Table 3. Table 3: Grid convergence analysis based on averaged Nusselt number on the hot wall.

No. of grid nodes	Averaged Nusselt number ($\overline{Nu}$)	Relative difference
$32\times 32$	11.3216	$28.2306\%$
$64\times 64$	9.2967	$5.2970\%$
$128\times 128$	8.9193	$1.0218\%$
$256\times 256$	8.84401	$0.1698\%$
$512\times 512$	8.82901	–

Table 4. Table 4: Effects of ε 𝜀 \varepsilon on L 𝐿 L and MAREs in mean and variance predictions for N s = 400 subscript 𝑁 𝑠 400 N_{s}=400 in the stochastic heat-driven cavity flow problem.

Tolerance		Size of selected POD basis		MARE in mean		MARE in variance
$\epsilon$		$L$		${\epsilon }_{r,\mu }^{FPCE}$		${\epsilon }_{r,var}^{FPCE}$
${10}^{-4}$		3		8.5751E-05		1.6665
${10}^{-5}$		9		9.0676E-05		2.4500E-2
${10}^{-6}$		12		9.0842E-05		1.3300E-02
${10}^{-7}$		14		3.0573E-05		5.6000E-03
${10}^{-8}$		18		2.5137E-05		5.4000E-03
${10}^{-9}$		28		2.3924E-05		5.4000E-03

Table 5. Table 5: Effects of N s subscript 𝑁 𝑠 N_{s} on L 𝐿 L and MAREs in mean and variance predictions for ε = 10 − 8 𝜀 superscript 10 8 \varepsilon=10^{-8} for the stochastic heat-driven cavity flow problem.

No. of snapshots		Size of selected POD basis		MARE in mean		MARE in variance
$N_s$		$L$		${\epsilon }_{r,\mu }^{FPCE}$		${\epsilon }_{r,var}^{FPCE}$
50		16		2.2670E-05		7.9990E-01
100		17		2.8925E-05		3.9400E-02
200		18		1.5711E-05		7.4000E-03
400		18		2.2587E-06		5.4000E-03
600		18		4.1937E-06		5.2000E-03

Table 6. Table 6: MAREs in mean and variance predictions for ε = 10 − 8 𝜀 superscript 10 8 \varepsilon=10^{-8} for different number of snapshots compared with the case of N s = 1000 subscript 𝑁 𝑠 1000 N_{s}=1000 .

No. of snapshots		Size of selected POD basis		MARE in mean		MARE in variance
$N_s$		$L$		${\epsilon }_{r,\mu }^{1000}$		${\epsilon }_{r,var}^{1000}$
200		18		3.1563E-05		2.4200E-02
400		18		7.1596E-06		6.1000E-03
600		18		3.8645E-06		1.5000E-03
800		18		2.6028E-06		7.9285E-04
1000		18		–		–

Table 7. Table 7: Grid convergence analysis based on the average value of temperature field on different meshes.

No. of grid nodes ($N$)	Average value of $T\left(𝒙\right)$ ($\frac{1}{N}\sum T(x,y)$)	Relative difference
$9\times 9$	0.2344	$6.16\%$
$17\times 17$	0.2461	$1.47\%$
$33\times 33$	0.2490	$0.29\%$
$65\times 65$	0.2498	–

Table 8. Table 8: Comparison of computational cost and performance of ROM-II and POD–PCE for the stochastic Ackley function.

Method	Size of DoE		RRMSE					Computational cost [s] (Excluding simulation cost)
			Minimum	Maximum	Mean	Standard deviation
ROM-II	400		3.2987E-04	3.0000E-03	8.1723E-04	2.6191E-04		2.73
ROM-II	10000		2.5521E-05	5.2281E-04	8.6826E-05	5.0008E-05		430
POD–PCE	400		5.9852E-06	5.7499E-04	7.7981E-05	5.8665E-05		5.27

Table 9. Table 9: Comparison of computational cost and performance of ROM-II and POD–PCE for the stochastic heat-driven cavity flow.

Method	Size of DoE		RRMSE					Computational cost [s] (Excluding simulation cost)
			Minimum	Maximum	Mean	Standard deviation
ROM-II	400		2.1777E-05	1.0900E-02	3.4103E-04	6.5481E-04		0.573
ROM-II	1000		1.0370E-05	1.0700E-02	2.6235E-04	5.4645E-04		1.78
POD–PCE	400		1.6801E-06	1.0200E-02	2.4889E-04	5.2770E-04		4.90

Table 10. Table 10: Comparison of computational cost and performance of ROM-II and POD–PCE for the stochastic heat diffusion.

Method	Size of DoE		RRMSE					Computational cost [s] (Excluding simulation cost)
			Minimum	Maximum	Mean	Standard deviation
ROM-II	400		0.0000E-00	1.3000E-03	4.0821E-04	2.6473E-04		0.161
ROM-II	2800		0.0000E-00	9.7399E-04	1.5943E-04	1.8109E-05		11.9
POD–PCE	400		0.0000E-00	9.7555E-04	1.5968E-04	1.8535E-05		4.93

Equations

$$\mathbf{U}={\left[U^{(1)},U^{(2)},\dots ,U^{(N)}\right]}^T,$$

$${\mathbf{U}}_{n_s}={\left[U_{n_s}^{(1)},\dots ,U_{n_s}^{(N)}\right]}^T,U_{n_s}^{(n)}=U({\mathbf{x}}^{(n)};{\mathbf{\xi }}^{(n_s)}),n=1,\dots ,N,n_s=1,\dots ,N_s.$$

$$U(\mathbf{x};\mathbf{\xi })\approx \sum _{l=1}^L{u}_l(\mathbf{\xi }){\psi }_l(\mathbf{x}),$$

$${\mathcal{UU}}^T\mathbf{\psi }=\lambda \mathbf{\psi },$$

$$\mathcal{U}={\mathbf{\Psi \Sigma }\mathbf{V}}^T,\mathbf{\Sigma }=\left(\begin{array}{cc}\mathbf{D}& 0\\ 0& 0\end{array}\right),$$

$$\mathcal{U}{\mathbf{\nu }}_{n_r}={\sigma }_{n_r}{\mathbf{\psi }}_{n_r}\text{and}{\mathcal{U}}^T{\mathbf{\psi }}_{n_r}={\sigma }_{n_r}{\mathbf{\nu }}_{n_r}\text{for}{n}_r=1,\dots ,N_r.$$

$${\mathcal{U}}^T\mathcal{U}{\mathbf{\nu }}_{n_r}={\sigma }_{n_r}^2{\mathbf{\nu }}_{n_r}\text{and}\mathcal{U}{\mathcal{U}}^T{\mathbf{\psi }}_{n_r}={\sigma }_{n_r}^2{\mathbf{\psi }}_{n_r}\text{for}{n}_r=1,\dots ,N_r.$$

$${\mathcal{U}}^T\mathcal{U}\mathbf{\nu }=\lambda \mathbf{\nu }.$$

$$\Vert \mathcal{U}-\tilde{\mathbf{\Psi }}{\tilde{\mathbf{\Psi }}}^T\mathcal{U}{\Vert }_F^2$$

$$=\underset{{\mathbf{\Psi }}^L\in {\mathbb{R}}^{N\times L}}{\min }\sum _{n_s=1}^{N_s}\Vert {\mathbf{U}}_{n_s}-{\mathbf{\Psi }}^L({\mathbf{\Psi }}^L{)}^T{\mathbf{U}}_{n_s}{\Vert }_2^2=\sum _{n_s=L+1}^{N_s}{\lambda }_{n_s}=\sum _{n_s=L+1}^{N_r}{\lambda }_{n_s},$$

$$\Vert \mathbf{A}{\Vert }_F=\sqrt{\sum _{n_1=1}^{N_1}\sum _{n_2=1}^{N_2}|A_{n_1,n_2}{|}^2}\text{for}\mathbf{A}\in {\mathbb{R}}^{N_1\times N_2}.$$

$$\frac{{\sum }_{l=1}^L{\lambda }_l}{{\sum }_{l=1}^{N_r}{\lambda }_l}>1-\epsilon ,$$

$$V_N=\left\{U(\mathbf{x};\mathbf{\xi }):\mathbf{x}\in \Omega ,\mathbf{\xi }\in \Gamma \right\},V_{N_s}=span\left\{{\mathbf{U}}_1,{\mathbf{U}}_2,\dots ,{\mathbf{U}}_{N_s}\right\},\Xi =span\left\{{\mathbf{\xi }}^{(1)},{\mathbf{\xi }}^{(2)},\dots ,{\mathbf{\xi }}^{(N_s)}\right\},$$

$$\mathbb{E}\left[{\varphi }_i^{(k)}({\xi }_k){\varphi }_j^{(k)}({\xi }_k)\right]={\int }_{{\Gamma }_k}{\varphi }_i^{(k)}({\xi }_k){\varphi }_j^{(k)}({\xi }_k){\rho }_k({\xi }_k)d{\xi }_k={\delta }_{ij},$$

$$u_l(\mathbf{\xi })=\sum _{\mathbf{\alpha }\in {\mathbb{N}}^d}{c}_{\mathbf{\alpha }}^{(l)}{\Phi }_{\mathbf{\alpha }}(\mathbf{\xi }),$$

$$u_l(\mathbf{\xi })\approx \sum _{\mathbf{\alpha }\in \mathcal{A}}{c}_{\mathbf{\alpha }}^{(l)}{\Phi }_{\mathbf{\alpha }}(\mathbf{\xi }),$$

$$\mathcal{A}=\left\{\mathbf{\alpha }\in {\mathbb{N}}^d:|\mathbf{\alpha }|\le p\right\},card\{\mathcal{A}\}=M=\frac{(p+d)!}{p!d!}.$$

$$\mathbb{E}\left[{\Phi }_{\mathbf{\alpha }}(\mathbf{\xi }){\Phi }_{\mathbf{\beta }}(\mathbf{\xi })\right]={\delta }_{\mathbf{\alpha }\mathbf{\beta }},$$

$$\mu (u_l)=\mathbb{E}\left[u_l(\mathbf{\xi })\right]\approx {\int }_{\Gamma }\left(\sum _{\mathbf{\alpha }\in \mathcal{A}}{c}_{\mathbf{\alpha }}^{(l)}{\Phi }_{\mathbf{\alpha }}(\mathbf{\xi })\right)\mathbf{\rho }(\mathbf{\xi })d\mathbf{\xi }=c_0^{(l)},$$

$$var(u_l)=\mathbb{E}\left[(u_l(\mathbf{\xi })-\mu (u_l){)}^2\right]\approx \sum _{0<|\mathbf{\alpha }|\le p}(c_{\mathbf{\alpha }}^{(l)}{)}^2.$$

$$u_l(\mathbf{\xi })=\sum _{m=0}^{M-1}{c}_k^{(l)}{\Phi }_m(\mathbf{\xi })+{\epsilon }_M=\mathbf{\Phi }{\mathbf{c}}^{(l)}+{\epsilon }_M,$$

$${\hat{\mathbf{c}}}^{(l)}=\arg \underset{{\mathbf{c}}^{(l)}}{\min }\mathbb{E}\left[{\left(\mathbf{\Phi }{\mathbf{c}}^{(l)}-u_l(\mathbf{\xi })\right)}^2\right].$$

$$\mathbb{E}\left[{\left(\mathbf{\Phi }{\mathbf{c}}^{(l)}-u_l(\mathbf{\xi })\right)}^2\right]$$

$$\approx \frac{1}{N_{\xi }}\sum _{n_{\xi }=1}^{N_{\xi }}{\left(\mathbf{\Phi }({\mathbf{\xi }}^{(n_{\xi })}){\mathbf{c}}^{(l)}-u_l({\mathbf{\xi }}^{(n_{\xi })})\right)}^2=\frac{1}{N_{\xi }}\Vert \mathbf{A}{\mathbf{c}}^{(l)}-{\mathbf{u}}_l{\Vert }_2^2.$$

$$A_{n_{\xi },m}={\Phi }_m({\mathbf{\xi }}^{(n_{\xi })}),n_{\xi }=1,\dots ,N_{\xi },\text{and}m=0,\dots ,M-1.$$

$${\hat{\mathbf{c}}}^{(l)}=({\mathbf{A}}^T\mathbf{A}{)}^{-1}{\mathbf{A}}^T{\mathbf{u}}_l.$$

$${\hat{\mathbf{c}}}^{(l)}=\arg \underset{{\mathbf{c}}^{(l)}}{\min }\mathbb{E}\left[{\left(\mathbf{\Phi }{\mathbf{c}}^{(l)}-u_l(\mathbf{\xi })\right)}^2\right]+\gamma \Vert {\mathbf{c}}^{(l)}{\Vert }_1,$$

$$U(\mathbf{x};\mathbf{\xi })\approx \tilde{U}(\mathbf{x};\mathbf{\xi })=\sum _{l=1}^L\left(\sum _{m=0}^{M-1}{c}_m^{(l)}{\Phi }_m(\mathbf{\xi })\right){\psi }_l(\mathbf{x}).$$

$$\mu (U(\mathbf{x};\mathbf{\xi }))=\mathbb{E}\left[U(\mathbf{x};\mathbf{\xi })\right]\approx \mathbb{E}\left[\sum _{l=1}^L{u}_l(\mathbf{\xi }){\psi }_l(\mathbf{x})\right]=\sum _{l=1}^L\mathbb{E}\left[u_l(\mathbf{\xi })\right]{\psi }_l(\mathbf{x})=\sum _{l=1}^L{c}_0^{(l)}{\psi }_l(\mathbf{x}),$$

$$var(U(\mathbf{x};\mathbf{\xi }))$$

$$Cov(u_{l_1},u_{l_2})$$

Full text

1 Introduction

Uncertainty quantification (UQ) is the process of quantifying the effects of input uncertainties on system responses. Recently, UQ has been widely explored in many fields of science and engineering, such as fluid [1, 2, 3] and structural mechanics [4, 5]. To reproduce the fine-scale structures of such physical phenomena, a large number of degree of freedom (DOF) for simulations is required, which clearly commensurates with increasing computational cost. Hence, for UQ, problems are in general too large to be tackled using standard techniques. This is because approximating system responses with a large number of random parameters on a fine mesh is usually computationally challenging. This emphasizes the need for establishing a model reduction method to reduce computational costs and storage requirements in UQ analysis of the problems with high-dimensional physical and random spaces.

A fundamental problem in UQ is approximating a computational model with random parameters. To this end, many numerical methods have been well developed in recent years (see [6, 7, 8, 9, 10, 11, 12] and references therein). One of the most widely used techniques in UQ is polynomial chaos expansion (PCE), which was first proposed by Ghanem and Spanos as finite-dimensional Wiener polynomial chaos [6, 13]. Later, Xiu and Karniadakis extended it to other random variables with basis functions from the Askey family of hypergeometric polynomials [14], which is known as generalized PCE. Furthermore, PCE was extended for random variables with arbitrary distributions in [15, 16]. To obtain deterministic coefficients of the PCE, two kinds of approaches are used: intrusive and non-intrusive methods. Intrusive methods need solvers for the resulting coupled deterministic equations, which can be very complicated if the form of the equations is nontrivial and nonlinear. However, in the non-intrusive methods, existing solvers are considered as ’black box’ without any modifications of the underlying computer codes. The PCE methods always exhibit high convergence rates with increasing expansion order, provided that the solutions are sufficiently smooth with respect to the random variables. However, a well-known shortcoming of these methods is the ”curse of dimensionality,” i.e., the simulation cost grows rapidly with increasing dimensions of random space. To mitigate this issue, various methods have been developed based on sampling strategies [17, 18, 19] and basis reduction [9, 12]. Nevertheless, for problems with high-dimensional physical space, constructing a PCE is still computationally challenging. This is because the computational cost depends not only on the cost of the PCE construction, but also on the number of spatial DOFs. An efficient reduced-order modeling (ROM) method is thus essential for such problems.

ROM [20] aims to approximate the original model using an accurate reduced-order model with a much smaller number of DOFs. Similar to the PCE methods, ROM methods can be classified into two categories in terms of dependency on the governing equations, namely, intrusive and non-intrusive [21]. Over the previous decades, many methods have been developed for ROMs (see [22, 23, 24] and references therein). In the present study, we focus on non-intrusive ROM methods based on proper orthogonal decomposition (POD) [26]. Usually, POD-based ROM methods construct a set of basis functions from the snapshots of the original full order model and then approximate the model responses with these basis functions. To this end, a variety of methods have been introduced for POD-based ROM. Audouze et al. [25, 27] proposed a non-intrusive POD-based ROM using radial basis functions (RBFs) for approximating solutions of nonlinear time-dependent parameterized partial differential equations. Xiao et al. [28] proposed a non-intrusive ROM based on constrained POD and the Kriging interpolation method, where they used Kriging interpolation to compute the POD coefficients. Guo and Hesthaven [29, 30] approximated the POD coefficients using Gaussian process regression in ROMs for nonlinear structural analysis and time-dependent problems. Noack et al. [31] introduced neural networks for ROMs. Later, a non-intrusive ROM method based on POD and RBF artificial neural networks (ANNs) was proposed by Vasile et al. [32]. Furthermore, ANNs, particularly multi-layer perceptrons, were introduced for ROMs to accurately approximate coefficients of the reduced model [33, 34]. Recently, a more extensive and comprehensive discussion on the use of machine learning for the approximation of POD coefficients was proposed in [35]. For more reviews of non-intrusive POD based ROM methods, please refer to [24]. In the present study, we focus on developing a non-intrusive ROM method for stochastic representations in UQ problems based on POD and PCE methods.

Recently, some algorithms have been proposed based on the POD and PCE methods (see [36, 37, 38, 39]). In [36], Doostan et al. proposed an intrusive stochastic model reduction method for creating an alternative optional basis for the PCE by adapting the Karhunen–Lovève (KL) expansion [40]. More precisely, a stochastic basis was first obtained on a coarse mesh by means of a KL expansion of the coarse-scale response and then used as the basis functions in fine discretization analysis. According to the results of the numerical examples in [4] and [36], dramatic reduction of computational cost was achieved with little loss in accuracy. However, the method is intrusive, which may be very complicated or even impossible to implement in some cases. To mitigate this problem, Raise et al. [37, 38] extended this idea to a non-intrusive PCE method, which is more efficient and broadly applicable. Moreover, to further reduce the computational cost of PCE construction, Abraham et al. [39] further extended this non-intrusive method by introducing a sparse polynomial chaos expansion method [12]. Clearly, the new non-intrusive sparse PCE-based model reduction scheme brings a significant speed-up for a wide range of engineering problems, especially for those with high-dimensional random space. However, a full PCE still needs to be constructed on a coarse discretization, which is a large computational challenge for some problems with high resolution in spatial discretization.

In this study, a non-intrusive ROM method is proposed for UQ problems with high-dimensional physical and random spaces. The method combines the advantages of both the POD and PCE methods. However, in comparison with [41], where POD was first constructed from a data set describing the ice shapes as a reduced-order model and then PCE was constructed based on this model to study the effects of uncertain parameters in POD representation, in our method, the POD basis is directly constructed for fine discretization using the method of snapshots [20]. Then, a sparse PCE method proposed in [9] is used to approximate the random coefficients of the POD representation. The method is referred to as POD–PCE. To investigate the performance of this method, three numerical examples were tested, namely, a highly irregular Ackley function with three random parameters, two-dimensional (2D) heat-driven cavity flow with a stochastic boundary temperature, and 2D heat diffusion with stochastic conductivity.

The rest of this paper is organized as follows. In Section 2, we present the POD–PCE method for the representation of the stochastic field, where the construction of the POD basis based on the method of snapshots is discussed, and a sparse PCE method is introduced to approximate the coefficients of the POD representation. Section 3 is dedicated to the numerical tests of the POD–PCE method. Finally, in Section 4 we provide our conclusions.

2 Representation of the random field

We consider a random-field model response $U(\bm{x};\bm{\xi})$, where $\bm{\xi}=(\xi_{1},\ldots,\xi_{d})\in\Gamma=\prod_{k=1}^{d}\Gamma_{k}\subset\mathbb{R}^{d}$ is a random vector with independent components described by the joint probability density function (PDF) $\bm{\rho}(\bm{\xi})=\prod_{k=1}^{d}\rho_{k}(\xi_{k})$, where $\rho_{k}(\xi_{k})$ is the PDF of the $k$-th component of the random vector. $\bm{x}$ is a spatial vector, which varies in a compact set $\Omega$ of $\mathbb{R}^{s}\ (s=1\ \textrm{or}\ 2\ \textrm{or}\ 3)$. In the present study, we focus on the case of $s=2$. Thus, the variable $\bm{x}$ is deterministic and varies on a grid of size $N$ corresponding to a discretization of $\Omega$. Therefore, $U(\bm{x};\bm{\xi})$ is an $N$-dimensional output given by

[TABLE]

where $U^{(n)}=U(\bm{x}^{(n)};\bm{\xi})$ and $n=1,2,\ldots,N$. Our objective is to develop a reduced-order model for $U(\bm{x};\bm{\xi})$. The procedure of the construction of the reduced-order model is described in the remainder of this section.

2.1 Proper orthogonal decomposition

We first consider a POD method to build a reduced-order model. POD is a powerful and elegant method of data analysis aimed at deriving low-order models of dynamic systems. It was first proposed by Kosambi [26], and was widely used in various fields of engineering and science under a variety of different names such as Karhunen–Loève decomposition [40] in stochastic process modeling and principal component analysis in statistical analysis [42, 43]. The POD method essentially provides an optimally ordered orthonormal basis in the least-squares sense for a given set of theoretical, experimental, or computational data of high-dimensional systems [44]. The optimal lower-dimensional approximations for the given data are then obtained by truncating the optimal basis. In our study, the POD basis vectors are computed empirically using sampled data collected over a subset of the random space $\Omega_{r}$, typically using the method of snapshots [20].

Consider an ensemble of snapshots, $\left\{\bm{U}_{1},\bm{U}_{2},\ldots,\bm{U}_{N_{s}}\right\}$, where $\bm{U}_{n_{s}}\in\mathbb{R}^{N}$ denotes the $n_{s}$-th snapshot, and is given by

[TABLE]

Let $\bm{\mathcal{U}}=\left[\bm{U}_{1},\bm{U}_{2},\ldots,\bm{U}_{N_{s}}\right]$ be a real-valued snapshot matrix of rank $N_{r}\leq min\{N,N_{s}\}$. Then, the POD representation of $U(\bm{x};\bm{\xi})$ is defined as

[TABLE]

where $\psi_{l}(\bm{x})$ are the most characteristic structure functions extracted from the spatial domain $\Omega$ based on the snapshot matrix $\bm{\mathcal{U}}$, $u_{l}(\bm{\xi})$ are the corresponding coefficients for a given random vector $\bm{\xi}$, and $L\in\{1,\ldots,N_{r}\}$ is the number of selected structure functions. To determine a proper orthonormal basis $\left\{\bm{\psi}_{l}\right\}_{l=1}^{L}\ (\bm{\psi}_{l}=\left[\psi_{l}(\bm{x}^{(1)}),\ldots,\psi_{l}(\bm{x}^{(N)})\right]^{T})$ of rank $L$ for the given snapshots, we need to solve the following eigenvalue problem

[TABLE]

which is strongly related to the singular value decomposition (SVD) [48] of matrix $\bm{\bm{\mathcal{U}}}$; i.e.,

[TABLE]

where $\bm{D}=diagonal(\sigma_{1},\ldots,\sigma_{N_{r}})\in\mathbb{R}^{N_{r}\times N_{r}}$ containing the singular values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{N_{r}}\geq 0$ and zeros denote matrices of appropriate dimensions with zero-elements. $\bm{\Psi}\in\mathbb{R}^{N\times N}$ and $\bm{V}\in\mathbb{R}^{N_{s}\times N_{s}}$ are unitary matrices, i.e., $\bm{\Psi}\bm{\Psi}^{T}=\bm{\Psi}^{T}\bm{\Psi}=\bm{I}_{N}$ and $\bm{VV}^{T}=\bm{V}^{T}\bm{V}=\bm{I}_{N_{s}}$. Let $\bm{\psi}_{n_{r}}\in\mathbb{R}^{N}$ and $\bm{\nu}_{n_{r}}\in\mathbb{R}^{N_{s}}$ be the $n_{r}$-th column vectors of $\bm{\Psi}$ and $\bm{V}$, respectively. Then we have

[TABLE]

Thus,

[TABLE]

Hence, $\left\{\bm{\psi}_{n_{r}}\right\}_{{n_{r}}=1}^{N_{r}}$ and $\left\{\bm{\nu}_{n_{r}}\right\}_{n_{r}=1}^{N_{r}}$ are eigenvectors of $\bm{\mathcal{U}}\bm{\mathcal{U}}^{T}$ and $\bm{\mathcal{U}}^{T}\bm{\mathcal{U}}$, respectively, with eigenvalues $\lambda_{n_{r}}=\sigma_{n_{r}}^{2}>0$, ${n_{r}}=1,\ldots,N_{r}$. Under the assumption of $N_{s}\ll N$, we reformulate the computation of POD modes as an $N_{s}\times N_{s}$ eigenvalue problem that is computationally efficient; i.e.,

[TABLE]

According to Eq. (6), $\left\{\bm{\psi}_{l}\right\}_{l=1}^{L}\ (L\in\{1,\ldots,N_{r}\},\ \bm{\psi}_{l}=\bm{\mathcal{U}}\bm{\nu}_{l}/\sqrt{\lambda_{l}})$ is the POD basis of rank $L$ that we need [47]. Then, a reduced-order model is obtained by projecting the full order model onto the space spanned by the POD basis of rank $L$. The POD basis is optimal in the sense that it minimizes the least-squares error of the snapshot representation for the reduced basis of size $L$. Let $\widetilde{\bm{\Psi}}\in\mathbb{R}^{N\times L}$ be a matrix with a pairwise orthonormal vector $\bm{\psi}_{l}$. Then, the projection error of the snapshots can be evaluated as

[TABLE]

where $\bm{\Psi}^{L}=\left[\bm{\psi}_{l_{1}},\ldots,\bm{\psi}_{l_{L}}\right]\ (l_{i},\in\{1,\ldots,N_{r}\},i=1,\ldots,L)$, $\|\cdot\|_{2}$ is the general $L_{2}$ norm for a vector (e.g.,$\|\bm{\psi}_{l}\|_{2}=\sqrt{\sum_{n=1}^{N}\psi_{l}^{2}(\bm{x}^{(n)})}$), and $\|\cdot\|_{F}$ denotes the Frobenius norm given by

[TABLE]

As seen in Eq. (9), the square of the error in the snapshot representation is defined by the sum of the eigenvalues of the eigenvalue problem (8) corresponding to the modes not included in the selected basis. Furthermore, a fast decay of the eigenvalues of $\bm{\mathcal{U}}^{T}\bm{\mathcal{U}}$ is usually assumed in POD. Therefore, the size of the POD basis can be chosen based on the size of the dominant eigenspace such that

[TABLE]

where $\varepsilon$ is a user-specified tolerance, often considered to be $10^{-3}$ or smaller [23]. The numerator of Eq. (11) is often referred to as the ”energy” captured by the POD modes.

Let

[TABLE]

where $V_{N}$ is the realization domain of $U(\bm{x;\xi})$, $V_{N_{s}}\subset V_{N}$ is spanned by the snapshots and $\Xi$ is a subset of $\Gamma$ spanned by a set of selected random vectors. Then, if the random vectors are selected sufficiently well, $U(\bm{x;\xi})$ can be approximated using the basis of $V_{N_{s}}$. As mentioned above, the main result of POD is that the optimal subspace $V_{L}$ of dimension $L\ll\min\{N,N_{s}\}$ representing the snapshots is given by $V_{L}=span\left\{\bm{\psi}_{1},\bm{\psi}_{2},\ldots,\bm{\psi}_{L}\right\}$. Thus, for any random vector $\bm{\xi}\in\Gamma$, $U(\bm{x;\xi})$ can be approximated as shown in Eq. (3). Clearly, the coefficients $u_{l}(\bm{\xi})$ in Eq. (3) vary with the random vector $\bm{\xi}$, which means that we can construct a mapping between the POD coefficients and the random vector $\bm{\xi}$. To this end, a sparse PCE method is introduced in the following subsection.

Remark 1**.**

The POD bases optimally approximate ensemble snapshots of field $U(\bm{x;\xi})$ in the least-squares sense, but they are not bases approximating the field $U(\bm{x;\xi})$. Only when the random vectors are selected well enough can the POD bases be used to approximate the field data. Hence, choosing good random vectors for the snapshots is important in the method of snapshots. In the present study, the snapshots are chosen by a sequential sampling scheme, as described in Section 2.4.

2.2 Polynomial chaos expansion

We now present the basic formulation of the stochastic representations of the coefficients $u_{l}(l=1,\ldots,L)$ in the POD representation of the random field $U(\bm{x};\bm{\xi})$ using the PCE method.

The basic idea behind PCE is to use orthogonal polynomials in terms of random variables to approximate functions of those random variables. As shown in [8, 14], the polynomial basis $\bm{\phi}^{(k)}(\xi_{k})(k=1,\ldots,d)$ can be chosen on the basis of the correspondence between the weighting function of the polynomial family and the PDF of the random variable $\xi_{k}$; for example, normal and uniform distributions correspond to Hermite and Legendre polynomial bases, respectively. Furthermore, the univariate polynomial bases $\bm{\phi}^{(k)}(\xi_{k})$ satisfy the following orthogonality:

[TABLE]

where $\delta_{ij}$ is the Kronecker delta function. Then, the PCE of $u_{l}(\bm{\xi})$ is defined as

[TABLE]

where $\Phi_{\bm{\alpha}}(\bm{\xi})=\phi_{\alpha_{1}}^{(1)}(\xi_{1})\cdots\phi_{\alpha_{d}}^{(d)}(\xi_{d})$ are multivariate polynomials orthonormal with respect to $\bm{\rho}(\bm{\xi})$ and $\bm{\alpha}\in\mathbb{N}^{d}$ is a multi-index that identifies the components of the multivariate polynomials $\Phi_{\bm{\alpha}}$ and satisfies that $|\bm{\alpha}|=\alpha_{1}+\cdots+\alpha_{d},\ |\bm{\alpha}|\geq 0$. In addition, $c_{\bm{\alpha}}^{(l)}$ are the deterministic coefficients corresponding to the polynomials. The convergence of the Eq. (13) should be justified for different polynomials in practice. In short, the Wiener–Hermite chaos [6] can be justified by the Cameron–Martin theorem [49]. Rigorous justifications about the convergence of generalized PCE can be found in [50].

Practically, the sum in Eq. (13) needs to be truncated to a finite sum by introducing the truncated polynomial chaos expansion:

[TABLE]

where $\mathcal{A}\subset\mathbb{N}^{d}$ is the set of selected multi-indices of multivariate polynomials which may have the following form:

[TABLE]

Because of the orthogonality relations in Eq. (12), it follows that the multivariate polynomials are also orthonormal:

[TABLE]

where $\delta_{\alpha\beta}$ is an extension Kronecker symbol to the multi-dimensional case. Based the orthogonality relations, we can easily estimate the statistical quantities of $u_{l}$; e.g., for a PCE with the first multivariate polynomial $\Phi_{\bm{0}}(\bm{\xi})=1,$ we have

[TABLE]

and

[TABLE]

To calculate the coefficients of the PCE for a given polynomial basis, many methods have been proposed, e.g., see [8, 9, 51, 52, 53]. In the present study, only regression-based non-intrusive methods are considered, i.e., the PCE coefficients are estimated by minimizing the mean squared error of the response approximation. A direct approach to estimating the coefficients is to set up a least-squares minimization problem [54]. Thus, Eq. (13) can be written as a sum of its truncated version and a residual:

[TABLE]

where $\varepsilon_{M}$ is the truncation error, $\bm{\Phi}=\left[\Phi_{0}(\bm{\xi}),\ldots,\Phi_{M-1}(\bm{\xi})\right]$ is the matrix that assembles the values of all the orthonormal polynomials in $\bm{\xi}$, and $\bm{c}^{(l)}=\left[c_{0}^{(l)},\ldots,c_{M-1}^{(l)}\right]^{T}$ is a vector containing the coefficients.

The least-squares minimization problem is then set up as:

[TABLE]

Given a sampling of size $N_{\xi}$ of the random vector $\{\bm{\xi}^{(n_{\xi})}\}_{n_{\xi}=1}^{N_{\xi}}$, we have $\bm{\mathcal{U}}=\left[\bm{U}_{1},\bm{U}_{2},\ldots,\bm{U}_{N_{\xi}}\right]$, where $N_{s}=N_{\xi}$. Based on Eq. (3), the corresponding coefficients can be estimated as $\bm{u}_{l}=\bm{\widetilde{\Psi}}^{T}\bm{\mathcal{U}}=\left[u_{l}(\bm{\xi}^{(1)}),\ldots,u_{l}(\bm{\xi}^{(N_{\xi})})\right]^{T},$ then we have

[TABLE]

where

[TABLE]

Then, the ordinary least-square solution of Eq. (20) is calculated as:

[TABLE]

According to Hosder et al. [52], $2M\sim 3M$ samples are sufficient for accurate approximation to the statistics at each polynomial degree. However, considering this design of experiment (DoE), the computational effort may become untenable for large values of $d$ or $p$; e.g., $d\geq 10$ or $p\geq 10$. To mitigate this issue, several methods have been developed (see [9, 12, 55] and references therein). In our study, the least-angle regression-based adaptive sparse PCE is adopted, which is developed based on the sparsity of effects; i.e., most computational models are driven by main effects and low-order interactions [9]. In fact, the sparse PCE solves a modified least-squares problem in Eq. (23) by adding a penalty term:

[TABLE]

where $\gamma$ is the regularization parameter and $\|\cdot\|_{1}$ is the $l_{1}$ norm that sums the absolute element of a vector.

Remark 2**.**

Conventional PCE methods focus on the problem with smooth model response with regard to input variables. Thus when there are steep gradients or finite discontinuities in the random space, these methods converge very slowly or even fail to converge due to a linear combination of polynomials. For these cases, local expansions may be more efficient than expansions with global polynomials, such as the wavelet expansions [56] and the multi-element generalized PCE [57, 58]. In the present study, we mainly focus on the problems with smooth responses with respect to input variables.

2.3 POD–PCE

Combining the POD representation in Eq. (3) and the PCE approximations in Eq. (14) of the undetermined coefficients of the POD basis, a whole representation of the random field $U(\bm{x};\bm{\xi})$ can be obtained as follows:

[TABLE]

This stochastic representation is referred to as POD–PCE in the present study. The method reduces the spatial dimensionality from $N$ to $L$. Furthermore, similar to the classical PCE methods, it can estimate the statistical quantities of $U(\bm{x};\bm{\xi})$ in terms of $\bm{\xi}$ using the orthogonality of the PCE; more precisely, the statistical quantities can be calculated by PCE coefficients and POD basis functions; e.g.,

[TABLE]

and

[TABLE]

where

[TABLE]

2.4 Design of experiment in POD–PCE

As presented in the above subsections, a good design of experiment (DoE) is necessary in both the POD and PCE parts. In the POD part, we have to generate a proper set of snapshots such that the space expanded by the POD basis extracted from these snapshots can approximate the random field response. In the PCE part, a proper number of random samples should be generated to obtain a convergent result under some accuracy requirements.

According to the discussion in [47], the continuous form of the eigenvalue problem (8) can be written as:

[TABLE]

where $\langle\cdot,\cdot\rangle_{\mathbb{R}^{N}}$ is the inner product in ${\mathbb{R}^{N}}$ and $\left\langle\bm{U},\bm{\psi}\right\rangle_{\mathbb{R}^{N}}=\bm{U}^{T}\bm{\psi}$. Then, we have

[TABLE]

Let

[TABLE]

We obtain the following eigenvalue problem:

[TABLE]

where $\bm{C}=\left(C_{n_{1},n_{2}}\right)_{N\times N}$ is the correlation matrix of the random field $U(\bm{x};\bm{\xi})$. In comparison with the eigenvalue problem (4), $\bm{\mathcal{U}\mathcal{U}}^{T}$ can be considered as an approximation to $\bm{C}$. The approximation accuracy depends on the method used to solve Eq. (32), such as the Latin hypercube sampling (LHS) method [7], sparse grid [59], and greedy sampling method [60]. In the present study, we approximate $\bm{C}$ by using the snapshots chosen by a sequential sampling scheme as follows:

Generate an initial set of random vectors as $\bm{\Theta}_{0}=\left\{\bm{\xi}^{(1)},\ldots,\bm{\xi}^{(N_{0})}\right\}$ using the LHS method [7]. 2. 2.

Calculate snapshots $\bm{U}_{n_{0}}$ ($n_{0}=1,\ldots,N_{0}$) and assemble the snapshot matrix $\bm{\mathcal{U}}=\left[\bm{U}_{1},\ldots,\bm{U}_{N_{0}}\right]$. 3. 3.

Solve the eigenvalue problem of $\bm{\mathcal{U}^{T}\mathcal{U}}$ to generate $\lambda_{n_{0}}^{0}$ and normalize the eigenvalues as $\hat{\lambda}_{n_{0}}^{0}=\frac{\lambda_{n_{0}}^{0}}{\lambda_{1}^{0}},\ \lambda_{1}^{0}\geq\cdots\geq\lambda_{N_{0}}^{0}$. 4. 4.

Enrich $\bm{\Theta}_{0}$ to be $\bm{\Theta}_{1}\supset\bm{\Theta}_{0}$ of size $N_{1}$ using a nested LHS method [45] and repeat (2)–(3) to obtain normalized eigenvalues $\left\{\hat{\lambda}_{n_{1}}^{1}\right\}_{n_{1}=1}^{N_{1}}$. 5. 5.

Select $N_{l}$ eigenvalues such that $\hat{\lambda}_{N_{l}}^{k}>\tau,\ k=0,1$ and compute the mean absolute relative error as $\varepsilon_{\lambda}=\frac{1}{N_{l}}\sum_{n_{l}=1}^{N_{l}}\left|\frac{\hat{\lambda}_{n_{l}}^{0}-\hat{\lambda}_{n_{l}}^{1}}{\hat{\lambda}_{n_{l}}^{1}}\right|$. 6. 6.

Update $\bm{\Theta}_{0}=\bm{\Theta}_{1},\ N_{0}=N_{1}$ and repeat (1)–(5) until $\varepsilon_{\lambda}<\hat{\tau}$. Then define the number of the snapshots as $N_{s}=N_{0}$ and select $\mathcal{D}_{POD}=\left\{\left(\bm{\xi}^{(n_{s})},\bm{U}_{n_{s}}\right),\ n_{s}=1,2,\dots,N_{s}\right\}$ as the DoE of the POD construction.

In the above procedure, tolerance $\tau$ is used to select the eigenvalues with large ”energy.” This is because the decay of the eigenvalue is usually rapid and the smallest eigenvalues have little influence on the total energy. The tolerance can be used to reduce the computation complexity of error $\varepsilon_{\lambda}$. For the tolerance $\hat{\tau}$, it is the stopping criterion of the procedure that ensures that the difference between the eigenvalues computed by the neighbor iterations is small enough. In the present study, we set $\tau=10^{-10}$ and $\hat{\tau}=0.05$ throughout the paper.

For the PCE part, the DoE is defined as $\mathcal{D}_{PCE}=\left\{\left(\bm{\xi}^{(n_{s})},\bm{\widetilde{\Psi}}^{T}\bm{\mathcal{U}}(\bm{\xi}^{(n_{s})})\right),\ n_{s}=1,2,\dots,N_{s}\right\}$ (see Section 2.2). This implies that we use the same size DoE for both POD and PCE. This is because we apply sparse PCE to approximate the coefficients of POD, which can obtain an accurate approximation by only a few samples of input variables. Additionally, the above snapshot selection procedure demonstrates that the final selected snapshots lead to a set of convergent normalized eigenvalues. This indicates that the eigenvalue has the same decay when there are $N_{s}$ or more snapshots. Thus, the $N_{s}$ snapshots can capture most of the information in the random field $U(\bm{x};\bm{\xi})$. Therefore, using the same DoE in PCE is reasonable. If the PCE approximation does not satisfy the accuracy requirement under the same DoE as that of POD, we can enrich the samples by using the nested LHS.

2.5 Error analysis of POD–PCE

For any given parameters $\bm{\xi}\in\Gamma$, we have

[TABLE]

Here, $\varepsilon_{L}(\bm{x})$ is the error caused by POD truncation and snapshot selection. Let $\bm{\varepsilon}_{M}=\left[\varepsilon_{M}^{(1)},\ldots,\varepsilon_{M}^{(L)}\right]^{T}$ and $\bm{\varepsilon}_{L}=\left[\varepsilon_{L}(\bm{x}^{(1)}),\ldots,\varepsilon_{L}(\bm{x}^{(N)})\right]^{T}$. Then, the error of POD–PCE approximation to $U(\bm{x;\xi})$ for a spatial discretization of $N$ grids is defined as follows:

[TABLE]

where $\widehat{\bm{\Psi}}\widehat{\bm{\Psi}}^{T}$ is the continuous covariance matrix of the random field, which has the form of Eq. (32), and $\|\cdot\|_{2,2}$ is a matrix norm induced by the $L_{2}$ norm, e.g.,

[TABLE]

Because $\widetilde{\bm{\Psi}}\widetilde{\bm{\Psi}}^{T}=\bm{I}_{N}$, $\left\|\widetilde{\bm{\Psi}}\right\|_{2,2}=1$. Then, we have

[TABLE]

where $\widetilde{C}$ is a constant related to the method for the snapshot selection, e.g., $\widetilde{C}\sim O\left(N_{s}^{-1/2}\right)$ for standard Monte Carlo methods [61] and $\widetilde{C}\sim O\left((N_{s})^{-k}(\log N_{s})^{(d-1)(k+1)}\right)$ for the Smolyak sparse-grid method [62], where $k$ is the degree of the polynomials. Therefore, $\varepsilon_{ss}$ mainly depends on the number of snapshots ($N_{s}$). As discussed in Section 2.1, $\varepsilon_{rb}$ depends on a user-specified tolerance $\varepsilon$; i.e., $\varepsilon_{rb}$ decreases as $\varepsilon$ decreases. $\varepsilon_{PCE}$ depends on the maximum order $p$ and the size of DoE ($N_{\xi}$) [8]. In the present study, the snapshots are chosen using a sequential sampling scheme, as shown in Section 2.4. The PCE construction uses a non-intrusive regression-based method, as proposed in [9]. Hence, the accuracy of POD–PCE depends on the number of the snapshots ($N_{s}$), tolerance $\varepsilon$, maximum polynomial order ($p$) and size of DoE ($N_{\xi}$).

2.6 Implementation of POD–PCE

In this section, we summarize the overall procedure followed to construct a POD–PCE representation for an $\mathbb{R}^{N}$-valued random field $U(\bm{x};\bm{\xi})$. We first choose a set of snapshots using a sequential sampling method as shown in Section 2.4. Second, we construct the POD basis using the selected snapshots and obtain a set of reduced basis in the following steps:

Assemble the snapshot matrix $\bm{\mathcal{U}}=\left[\bm{U}_{1},\ldots,\bm{U}_{N_{s}}\right]$ and solve the eigenvalue problem of $\bm{\mathcal{U}}^{T}\bm{\mathcal{U}}$ by Eq. (8) to generate $\lambda_{n_{s}}$ and $\bm{\nu}_{n_{s}}$, $n_{s}=1,2,\ldots,N_{s}$. 2. 2.

Set $\bm{\psi}_{n_{r}}=\frac{\bm{\mathcal{U}}\bm{\nu}_{n_{r}}}{\sqrt{\lambda_{n_{r}}}}$, $n_{r}=1,2,\ldots,N_{r}$. 3. 3.

Reduce the basis to be $\left\{\bm{\psi}_{l}\right\}_{l=1}^{L}$ such that condition (11) is satisfied. 4. 4.

Define the reduced POD basis matrix as $\widetilde{\bm{\Psi}}=\left[\bm{\psi}_{1},\bm{\psi}_{2},\ldots,\bm{\psi}_{L}\right]$.

Third, a sparse PCE method is utilized to approximate the undetermined coefficients of POD:

Prepare the DoE as $\left\{\left(\bm{\xi}^{(n_{s})},u_{l}(\bm{\xi}^{(n_{s})})\right),\ n_{s}=1,2,\dots,N_{s}\right\}$. 2. 2.

Select a set of orthogonal polynomials corresponding to the distribution of the random variables. 3. 3.

Construct the PCEs of $u_{l}(\bm{\xi})$ as Eq. (13) by solving the $l_{1}$ minimization problem (23).

Finally, a POD–PCE representation of $U(\bm{x};\bm{\xi})$ is constructed using Eq. (24).

All the PCE constructions in this paper are done using the least angle regression algorithm with UQLab [63].

3 Numerical examples

3.1 Highly irregular Ackley function

First, we consider a highly irregular Ackley function with three random parameters. The formulation of this function is given by

[TABLE]

where $\bm{\xi}=[\xi_{1},\xi_{2},\xi_{3}]$ is a random vector with independent components that are uniformly distributed over $[-1,1]^{3}$.

According to the discussion in [37], a discretization with $400\times 400$ nodes is required to reproduce the fine-scale structure of the deterministic Ackley function. Hence, such a fine mesh is adopted in the present study for UQ analysis of the stochastic Ackley function. Figure 1 shows the Ackley function on the mesh for a randomly generated parameter vector, where the structure of the function is well depicted.

As discussed in Section 2, to construct a POD–PCE representation, a proper set of snapshots of a size of $N_{s}=400$ was first selected by following the procedure described in Section 2.4. Figure 2 shows the normalized eigenvalues calculated using different numbers of snapshots in the snapshot selection procedure. We observe that there is a sharp drop in the first seven eigenvalues. This indicates that most of the information in $u(\bm{x;\xi})$ can be captured using the first seven modes. This can be confirmed by Table 1, where the result shows that the first seven modes capture $99.99\%$ of the total energy.

With the selected snapshots and corresponding random vectors, we constructed both POD–PCE and full PCE (FPCE) for the output field data of the stochastic Ackley function. Similar to the PCE construction for the POD coefficients in Section 2.2, the full PCE of $u(\bm{x;\xi})$ is constructed as follows:

[TABLE]

The maximum order of the polynomials in the PCE construction parts was set as $p=13$ according to the recommendation in [37]. Furthermore, Monte Carlo (MC) simulations were performed using $3\times 10^{3}$ LHS samples to validate the accuracy of POD–PCE. Given $N_{s}=400$, we investigated how accurately POD–PCE predicts the statistical quantities by changing the value of $\varepsilon$. The results are presented in Table 1, where the number of the selected POD bases (L) and the mean absolute relative errors (MAREs) in mean and variance predictions are presented for different values of $\varepsilon$. The MAREs between POD–PCE and full PCE are defined as follows:

[TABLE]

Similar definitions can be obtained for the MAREs between POD–PCE and MC; i.e., $\varepsilon_{r,\mu}^{MC}$ and $\varepsilon_{r,var}^{MC}$. The MAREs between POD–PCE and MC (or full PCE) in mean and variance predictions decrease and the number of the selected POD bases increases as the value of $\varepsilon$ decreases from $\varepsilon=10^{-4}$ to $\varepsilon=10^{-8}$. Furthermore, the decreasing feature tends to flat as $\varepsilon$ becomes smaller. In particular, when $\varepsilon\leq 10^{-8}$, the MAREs between POD–PCE and MC are nearly unchanged and those between POD–PCE and full PCE are also very small. This indicates that POD–PCE and PCE have almost the same accuracy in mean and variance predictions. However, it is worth noting that POD-PCE only needs to construct $40$ PCEs for the undetermined coefficients of the selected POD bases for $\varepsilon=10^{-8}$, while full PCE needs to construct $N=160000$ PCEs for all the field data. Clearly, an enormous amount computational cost is reduced using POD–PCE method. In summary, POD–PCE constructed based on $\varepsilon=10^{-8}$ can predict mean and variance with an accuracy similar to that of the classical full PCE economically.

Next, we considered the effects of the number of snapshots ($N_{s}$) on how accurately POD–PCE predicts the statistical quantities for the tolerance $\varepsilon=10^{-8}$. We are interested in investigating whether the POD–PCE constructed using a small number of snapshots can provide estimations of statistical quantities as reasonable as those obtained using the classical full PCE method. Table 2 shows the effects of $N_{s}$ on the size of the selected POD basis and MAREs in mean and variance predictions for the stochastic Ackley function. The results of the POD–PCE are agreement with those of full PCE as the number of snapshots increases. Furthermore, we observe that even the POD–PCE constructed using $50$ snapshots can predict the mean and variance with reasonable accuracy. However, the error in variance prediction is always larger than that in mean prediction, which implies that higher order moment prediction requires higher accuracy PCE; further, attention should be paid to the high order moment prediction for PCE construction.

Finally, the efficiency and accuracy of the method of snapshots is validated by the comparison of eigenvalues calculated using full PCE and the method of snapshots. In the method of snapshots, eigenvalues are calculated using $400$ snapshots by Eq. (8). For the full PCE method, the eigenvalues are obtained by solving the eigenvalue problem of the correlation matrix $\bm{R}=(r_{n_{1},n_{2}})_{N\times N}$ with

[TABLE]

where $c_{k}(\bm{x}^{(\cdot)})$ are deterministic coefficients of Eq. (39) corresponding to grid $\bm{x}^{(\cdot)}$. The comparison result is shown in Figure 3, where we only present eigenvalues corresponding to the modes truncated by the tolerance $\varepsilon=10^{-8}$. For the full PCE method, $42$ bases are selected, which is slightly larger than that in the method of snapshots. However, the comparison shows good agreement for almost all eigenvalues. Only tiny discrepancies exist in the calculation of the smallest eigenvalues. Therefore, the method of snapshots is accurate enough to be used to construct the POD basis compared with the methods discussed in [36, 37, 39].

3.2 Heat-driven cavity flow with a stochastic boundary temperature

To further investigate the performance of the proposed POD-PCE method, we consider a thermally driven cavity flow with stochastic boundary temperatures in a square domain, which is widely used as a numerical example to validate UQ algorithms (see [64, 65, 66] and references therein). Figure 4 shows the schematic of the problem with the computational domain, boundary conditions, and a realization of the stochastic temperature field. As described in [65, 66], the boundary temperature $T_{h}$ on the left vertical wall is uniform, while the boundary temperature $T_{c}<T_{h}$ on the right vertical wall is a spatially varying random field. In terms of dimensionless variables, the time-dependent governing equations for fluids are presented as follows:

[TABLE]

where ${\bf u}$, $p$, and $T$ are the velocity vector, pressure, and temperature, respectively. $\bm{e}_{y}$ denotes the unit normal vector in the vertical dimension. The dimensionless variables employed in the governing equations are defined as follows:

[TABLE]

[TABLE]

where $x$ (and $y$), $L=1$, $u_{i}$, $U_{0}=\sqrt{g\beta L(T_{h}-\overline{T}_{c})\alpha/\nu}$, $T_{h}$, and $\overline{T}_{c}$ respectively represent the Cartesian coordinates, characteristic length, velocity components, characteristic velocity, hot wall temperature, and cold wall mean temperature. In addition, $\alpha$, $\beta=0.5$, $\nu$, $\rho$, and $g$ respectively represent the thermal diffusivity, thermal expansion coefficient, kinematic viscosity, fluid density, and gravitational acceleration. The symbol  $\widetilde{}$  indicates a dimensional variable. Here, the Prandtl and Rayleigh numbers are set to $Pr=0.71$ and $Ra=10^{6}$, respectively.

The temperature on the hot wall is set to a constant value, i.e., $T_{h}=0.5$. On the cold wall, we consider a boundary temperature with stochastic fluctuations as follows:

[TABLE]

where $\overline{T}_{c}=-0.5$ is a constant mean temperature and $T_{c}^{{}^{\prime}}(y)$ is the random component, which can be expressed in a truncated Karhunen–Loève expansion as described in [65, 66]:

[TABLE]

where $\xi_{i}$ are uniform random variables assumed to be independently and uniformly distributed over $[-1,1]$ and $\lambda_{i}$ and $\psi_{i}(y)$ are the eigenvalues and eigenfunctions of the exponential correlation function

[TABLE]

with $\sigma_{T}=11/100$ and $l=1/21$. Here, $d=20$ means that only largest eigenvalues and the corresponding eigenfunctions are used to represent the random field $T_{c}^{{}^{\prime}}(y)$ with the random variables $\xi_{i}$.

By Eqs. (45) and (46), deterministic boundary temperatures on the cold wall can be obtained for any random vector $\bm{\xi}=[\xi_{1},\ldots,\xi_{d}]$. Then, computations are performed using the deterministic monolithic projection method developed in [67, 68, 69]. In this problem, a fine mesh is found by performing grid convergence analysis on the averaged Nusselt number $\overline{Nu}$ on the hot wall for meshes with $32\times 32$, $64\times 64$, $128\times 128$, $256\times 256$, and $513\times 513$ grid nodes. $\overline{Nu}$ is obtained by integrating the local Nusselt number $Nu=-(\partial T/\partial x)|_{x=0}$ along the hot wall. The results are present in Table 3. Finally, we chose the mesh with $256\times 256$ grid nodes as the fine mesh because it featured the smallest and most reasonable relative difference from the mesh with $513\times 513$ grid nodes. To investigate the performance of POD–PCE, the output temperature field denoted by $T(\bm{x;\xi})$ was considered in this numerical test.

Similar to the former example, the size of the snapshots was first chosen as $N_{s}=400$ using the sequential sampling scheme described in Section 2.4. Figure 5 presents all the normalized eigenvalues calculated by the considered snapshots in decreasing order. Compared with the stochastic Ackley function problem, the current calculated eigenvalues decay at a faster rate. Furthermore, we observed a sharp drop in the first three eigenvalues, which means that the first three modes capture most of the information in the temperature field.

Using $N_{s}=400$ snapshots, PCEs are constructed with $p=4$ for both POD–PCE and full PCE for $T(\bm{x;\xi})$. The maximum order of polynomials is chosen as that in [65], where a good performance of a sparse PCE with $p=4$ is shown in the representation of the stochastic solutions of Eqs. (42)–(44). Further, the fourth-order classical full PCE is constructed using $400$ snapshots and their corresponding random vectors throughout this example. The effects of $\varepsilon$ and $N_{s}$ on the accuracy of POD–PCE in the prediction of statistical quantities are presented in Tables 4 and 5, respectively. Similar to the previous example, there is a good agreement between POD–PCE results and full PCE results. In particular, the MAREs in mean predictions are very small, even for a large tolerance or fewer snapshots. This is because the first mode is the mean. Thus, if the first mode is selected, we can give a reasonable and accurate estimate of the mean value. For the MAREs in variance predictions, both $\varepsilon$ and $N_{s}$ have positive effects in that the MAREs decrease with decreasing $\varepsilon$ or increasing $N_{s}$, and finally, POD–PCE and full PCE obtain similar variance predictions. These indicate that POD–PCE has a similar accuracy to the classical full PCE in the predictions of statistical quantities of $T(\bm{x;\xi})$. Furthermore, it can reasonably predict the statistical quantities with fewer snapshots; e.g., using $N_{s}=100$ snapshots. In addition, one can find that the MARE of mean prediction increases when $N_{s}=600$. This is because larger input samples improve the accuracy of POD–PCE. However, the full PCE only used $400$ input samples; the MAREs therefore become larger than for when $N_{s}=400$. In fact, the mean prediction gradually converges as $N_{s}$ increases. This is illustrated in Table 6, where $\varepsilon_{r,\mu}^{1000}$ and $\varepsilon_{r,var}^{1000}$ denotes the MAREs computed for $N_{s}=1000$.

To further investigate the accuracy of POD-PCE, we first compared the prediction results of POD–PCE and full PCE for a realization of $T(\bm{x;\xi})$. The results are shown in Figure 6 with the result directly simulated from Eqs. (42)–(44). We observe that the contours of temperature fields are similar, and it is difficult to distinguish differences among these contour plots. This indicates that both full PCE and POD–PCE can give an accurate prediction of $T(\bm{x;\xi})$ and the errors in the full PCE and POD–PCE predictions are comparable. However, we guess that a few larger errors in the predictions of POD–PCE and full PCE may exist near the right cold wall due to the uncertain temperature boundary condition.

Next, we measured the error in the predictions of POD–PCE and full PCE by relative root mean-squared error (RRMSE), defined by

[TABLE]

where $HF$ denotes the high-fidelity result directly simulated from the full order model and $LF$ represents the low-fidelity result estimated from POD–PCE or full PCE. The results are shown in Figure 7, where we compare the distributions of RRMSE for the predictions of $T(\bm{x;\xi})$ in POD–PCE and full PCE. The RRMSEs are computed based on $100$ realizations of $T(\bm{x;\xi})$ calculated via direct simulation, POD–PCE, and full PCE for $100$ randomly generated input variable vectors $\{\bm{\xi}^{(k)}\}_{k=1}^{100}$. Large RRMSEs are mainly distributed near the right cold wall, especially at the right bottom angle of the computational domain. Furthermore, larger RRMSEs exist in the predictions of POD–PCE compared with those in full PCE. This is because of the error introduced by the truncated POD in POD–PCE method. In fact, these larger RRMSEs may tend to the values of the corresponding RRMSEs of the full PCE when we decrease the selection tolerance, which can be illustrated in Figure 8(b). In addition, we can find that the values of RRMSE of POD–PCE and full PCE are comparable in the predictions of $T(\bm{x;\xi})$. This indicates that the POD–PCE has an accuracy comparable to that of full PCE in the prediction of $T(\bm{x;\xi})$.

Finally, we investigate the effects of the number of snapshots ($N_{s}$) and the tolerance ($\varepsilon$) of POD basis selection on the accuracy of POD–PCE using the averaged RRMSEs over the temperature field. By Eq. (48), the averaged RRMSEs are given by

[TABLE]

For a fixed tolerance $\varepsilon=10^{-8}$, the averaged RRMSEs of $T(\bm{x;\xi})$ predications in the POD–PCE method are calculated for the aforementioned $100$ realizations for different $N_{s}$. The results are shown in Figure 8(a). The error decreases with the increasing number of snapshots. However, the error decreasing rate becomes smaller simultaneously. In general, when $100<N_{s}<400$, POD–PCE has a comparable averaged RRMSE to that of full PCE and when $N_{s}>400$, a similar averaged RRMSE to that of full PCE can be obtained. The difference between two average RRMSEs at $N_{s}=400$ is based on POD truncation. Therefore, the POD truncation tolerance ($\varepsilon$) should be decreased to improve the accuracy of POD–PCE. This is studied in Figure 8(b), where the averaged RRMSEs are computed for a fixed number of snapshots; i.e., $N_{s}=400$. Clearly, we observe that the error decreases as $\varepsilon$ decreases. Furthermore, we find that, when $\varepsilon\leq 10^{-8}$, POD–PCE can have similar averaged RRMSE to that of full PCE. In summary, if we use the same expansion order in the PCE construction parts for both POD–PCE and full PCE, the accuracy of POD–PCE depends on $N_{s}$ and $\varepsilon$. By accurately choosing the number of snapshots and tolerance value, the proposed POD–PCE method can deliver accuracy similar to that of the classical full PCE method.

3.3 Two-dimensional heat diffusion with stochastic conductivity

In this section, we consider a two-dimensional (2D) steady-state heat diffusion in a square plate. The governing equation is described as follows:

[TABLE]

where $k(x,y;\bm{\xi})$ is assumed to be a 2D random field defined as a function of the spatial vector $\bm{x}=(x,y)$ as well as the random vector $\bm{\xi}$. The schematic of the computational domain, boundary conditions, and a realization of the thermal conductivity are shown in Figure 9. The top boundary of the plate is at a high temperature $T_{h}$ and other boundaries of the plate are at a low temperature $T_{c}$. The thermal conductivity $k(\bm{x};\bm{\xi})$ is defined over $\mathcal{D}=[0,1]^{2}$ with known mean $\bar{k}(\bm{x})$ and covariance function:

[TABLE]

where $\sigma$ is standard deviation on the conductivity and $l_{1}$ and $l_{2}$ are the correlation lengths in $x$ and $y$ directions, respectively.

According to KL expansion, we can approximate $k(\bm{x;\xi})$ with $d$ inputs random variables as follows

[TABLE]

where $d$ is the truncation term and $\lambda_{n}$ and $\phi_{n}(\bm{x})$ are respectively eigenvalues and eigenfunctions obtained by solving the following 2D integral equation:

[TABLE]

For any given random vector $\bm{\xi}$, $k(\bm{x;\xi})$ can be approximated by Eq. (52). The deterministic computation of Eq. (50) is then solved using a centered finite difference scheme with Jacobi iterative method on a uniform grid. In the present study, the boundary temperatures are set to $T_{h}=1$ and $T_{c}=0$, the thermal conductivity mean and standard deviation are respectively set to $\bar{k}(\bm{x})=1$ and $\sigma=0.2$, and the correlation lengths are assumed to be $l_{1}=l_{2}=6$. In addition, the number of grid nodes is set to $N$ and the stopping criterion in the deterministic solver is set to be $10^{-10}$.

Table 7 shows the result of the grid convergence analysis based on the average values of temperature field on different meshes. It is clearly observed that the mesh with grid nodes $N=33\times 33$ is fine enough to be used to simulate the present heat diffusion problem. Thus, the following analyses are performed based on this mesh.

In this example, we focused on investigating the accuracy of the proposed POD–PCE method in the prediction of $T(\bm{x;\xi})$. Similar to the previous examples, $N_{s}=400$ snapshots was first chosen using the selection procedure in Section 2.4. The POD truncation criterion was defined as $\varepsilon=10^{-9}$ and $L=17$ POD bases were then selected based on this criterion. Next, we constructed the PCEs for both POD–PCE and full PCE methods using the DoE of the size $N_{s}=400$ as described in Section 2.6. Furthermore, all the PCEs are built using the third-order Legendre polynomials ($p=3$).

For a randomly generated thermal conductivity, we estimated $T(\bm{x;\xi})$ using POD–PCE and full PCE, respectively. The results are shown in Figure 10, compared with the results directly simulated using the full order model. In general, we observe a good agreement among these results. This indicates that the classical full PCE and POD–PCE both provide reasonable estimations for $T(\bm{x;\xi})$. To further investigate the prediction accuracy of the proposed POD–PCE method, we measured the error in the predictions of POD-PCE and full PCE by the RRMSE as defined in Eq. (48). Figure 11 shows the distributions of RRMSE for the predictions of $T(\bm{x;\xi})$ in full PCE and POD–PCE. The RRMSEs are calculated based on $100$ randomly generated thermal conductivity fields. Similar to heat-driven cavity flow, the distributions of the two RRMSE fields are similar, whereas more larger values are found in the RRMSE field for the POD–PCE method, especially at the bottom corners of the domain. This implies that two methods have similar accuracy. However, the performance of the full PCE is slightly better than that of the POD–PCE method. This is induced by the errors coming from the POD construction. We can improve the accuracy of the POD–PCE method by increasing the number of snapshots ($N_{s}$) or decreasing the POD truncation tolerance ($\varepsilon$). This is illustrated in Figure 12, where the averaged RRMSEs of $T(\bm{x;\xi})$ predictions using POD–PCE method for different number of snapshots and truncation tolerances are shown. The results in Figure 12(a) are calculated based on a fixed $\varepsilon=10^{-9}$, while the results in Figure 12(b) are calculated based on a fixed $N_{s}=400$. From the figure, we observe that the error is continuously decreasing for a fixed $\varepsilon=10^{-9}$ when $N_{s}>400$ though the decreasing rate is very small. In particular, for $N_{s}=400$, we find that the averaged RRMSE of POD–PCE is similar to that of full PCE when $\varepsilon=10^{-10}$. This re-emphasizes that POD–PCE can reach a similar accuracy as full PCE when we use the same size DoE in both methods.

Furthermore, we compared our method with the reduced-order modeling method (physical-base ROM-II) proposed in [25] from the view of computational cost using the above three examples. The results are presented in A. It is found that the physical-based ROM-II is usually faster than POD–PCE under the same DoE. However, the accuracy of the physical-based ROM-II is worse than that of the POD–PCE. To obtain a similar accuracy as the POD-PCE method, more simulations of the full order model are needed leading to higher computational cost compared with the POD–PCE method. In addition, the proposed POD–PCE method is much easier on the estimation of statistical quantities than the physical-base ROM-II method owing to its orthogonal property. Thus, the proposed POD–PCE is more efficient than the ROM-II method in the prediction of model outputs.

4 Conclusion

We developed a non-intrusive reduced-order modeling method for the stochastic representations in an uncertainty quantification analysis of problems with high-dimensional physical and random spaces. The method is based on a model reduction technique POD and an efficient stochastic spectral method PCE. POD is used to extract a set of reduced basis functions from a data ensemble collected from the full order model; i.e., snapshots. PCE is used to approximate the undetermined coefficients of the reduced basis. In particular, a non-intrusive sparse PCE is used to deal with highly random dimensionality problems. The proposed method inherits advantages of both POD and PCE methods. It can provide a reduced-order representation over the physical space as well as efficiently estimate the statistical quantities without increased computational cost. Three examples were tested to investigate the performance of the proposed method, i.e., an analytic highly irregular Ackley function with three random parameters, a heat-driven cavity flow with a stochastic boundary temperature, and a 2D heat diffusion with stochastic conductivity. The results show that the proposed method can reduce the model order from properly chosen snapshots at a reasonable accuracy. Furthermore, under the same parameter setting for the PCE construction, the proposed method can predict the model response and its statistical quantities with an accuracy similar to that of the classical full PCE method economically. In addition, the proposed method is more efficient compared with the existing radial basis function-based ROM-II proposed in [25]. These imply great potential for the proposed method applied in an uncertainty quantification analysis.

Acknowledgements

This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(Ministry of Science and ICT) (NRF-2017R1E1A1A0-3070161 and NRF-20151009350).

Appendix A Comparison between POD–PCE and ROM-II

This section provides a comparison between the proposed POD–PCE method and ROM-II method proposed in [25]. To this end, three examples in Section 3 are tested. The RRMSE defined in Eq. (48) is estimated as the measure of comparison. The results are shown in Tables 8–10 and Figures 13–15.

From the presented results, one can learn that under the same DoE of size $N_{s}=400$, the computational costs of the ROM-II for all the examples are much lower than the proposed POD–PCE, however, the accuracy performances are just the opposite. To obtain a comparable accuracy as the POD–PCE method, more simulations of the full order model are needed. For the Ackley function, $9600$ more simulations are needed, which considerably increases the computational cost, even for the ROM-II itself, the computational time excluding the simulation cost is much larger than that of the POD–PCE method. Similar observation can be made for the other two examples. For heat-driven cavity flow and heat-diffusion problem, $600$ and $2400$ more simulations are needed, respectively. Furthermore, one simulation cost of these two problems is much higher compared with the approaches cost. Therefore, the proposed POD–PCE is more efficient compared with the ROM-II method.