Gaussian Processes and Polynomial Chaos Expansion for Regression Problem: Linkage via the RKHS and Comparison via the KL Divergence

Abstract

In this paper, we examine two widely-used approaches, the polynomial chaos expansion (PCE) and Gaussian process (GP) regression, for the development of surrogate models. The theoretical differences between the PCE and GP approximations are discussed. A state-of-the-art PCE approach is constructed based on high precision quadrature points; however, the need for truncation may result in potential precision loss; the GP approach performs well on small datasets and allows a fine and precise trade-off between fitting the data and smoothing, but its overall performance depends largely on the training dataset. The reproducing kernel Hilbert space (RKHS) and Mercer’s theorem are introduced to form a linkage between the two methods. The theorem has proven that the two surrogates can be embedded in two isomorphic RKHS, by which we propose a novel method named Gaussian process on polynomial chaos basis (GPCB) that incorporates the PCE and GP. A theoretical comparison is made between the PCE and GPCB with the help of the Kullback–Leibler divergence. We present that the GPCB is as stable and accurate as the PCE method. Furthermore, the GPCB is a one-step Bayesian method that chooses the best subset of RKHS in which the true function should lie, while the PCE method requires an adaptive procedure. Simulations of 1D and 2D benchmark functions show that GPCB outperforms both the PCE and classical GP methods. In order to solve high dimensional problems, a random sample scheme with a constructive design (i.e., tensor product of quadrature points) is proposed to generate a valid training dataset for the GPCB method. This approach utilizes the nature of the high numerical accuracy underlying the quadrature points while ensuring the computational feasibility. Finally, the experimental results show that our sample strategy has a higher accuracy than classical experimental designs; meanwhile, it is suitable for solving high dimensional problems.

1. Introduction

Computer simulations are widely used in learning tasks, where a single simulation is an instance of the system [1,2]. A simple approach to the learning task is to randomly sample input variables and run the simulations for each input to obtain the features of the systems. Similar approaches are utilized in Monte Carlo techniques [3]. However, even a single simulation can be computationally costly due to its high complexity, and so, obtaining a trustworthy result via sufficient simulations becomes intractable. Mathematical methods and statistical theorems are introduced to generate surrogate models to replace the simulations, especially when dealing with complex systems with many parameters [4,5]. Although the main drawback of surrogate models is that only approximations can be obtained, they are computationally efficient whilst maintaining the essential information of the systems, hence analyzing the properties of the system. Attempting to construct surrogate models with an acceptable number of simulations necessitates the development of robust techniques to determine their reliability and validity [6,7,8]. Plenty of researchers are working on improving sampling strategies to decrease the number of simulations, which makes the task more significant [9]. With an increasing number of surrogate models being developed, there needs to be a comprehensive understanding of the uncertainties introduced by those models. The main purpose of uncertainty quantification (UQ) is to establish a relationship between input and output, i.e., the propagation of input uncertainties, and then to quantify the difference between surrogate models and original simulations. UQ can provide a measure of the surrogate model’s accuracy and an indication of how to update the model at the same time [10,11,12].

Denote f as a function (or simulator) of the original system, then given experimental design $\mathit{X}$, the output $Y=f\left(\mathit{X}\right)$ is produced, where caption notation is used because the input and output are usually vectors (or matrices) in simulations. From a statistical perspective, the input uncertainties are introduced by their randomness, so we represent the input with a random variable $\mathit{x}$, whose prior probability density function (PDF) is $p\left(\mathit{x}\right)$, such as the multivariate Gaussian distribution; as for the output uncertainties, a common technique is to integrate the system uncertainty and the approximation error as a noise term $\mathit{ϵ}$. In fact, the output $\mathit{y}$ is also a random variable $\mathit{y}=f\left(\mathit{x}\right)+\mathit{ϵ}$ determined by $f,\mathit{x}$ and $\mathit{ϵ}$. Now, suppose a surrogate $\overline{f}\left(\mathit{x}\right)$ is constructed to approximate $f\left(\mathit{x}\right)$, then UQ is used to identify the distribution and statistical features (for example, Kullback–Leibler divergence) of $\mathit{y}$, which are essential to the validation and verification of surrogates. Basically, there are two preconditions that need to be satisfied: firstly, the surrogate models are well defined, i.e., any $\overline{f}$ is a measurable function with respect to (w.r.t) corresponding probability space $p\left(\mathit{x}\right)$; secondly, techniques are needed that learn from the prior information to obtain the best guess of the true function.

There is a number of studies proposing different surrogates for specific applications in the literature, such as multivariate adaptive regression splines (MARS) [13], support vector regression (SVR) [14], artificial neural network (ANN) [15] for reliability and sensitivity analyses and kriging [16] for structural reliability analysis. We mainly focus on two popular methods that have been extensively studied recently. One popular method that is extensively studied in the literature is the polynomial chaos expansion (PCE), also known as a spectral approach [17]. PCE aims to represent an arbitrary random variable of interest as a spectral expansion function of other random variables with prior PDF. Xiu et al. [18,19,20] have generalized the PCE in terms of the Askey scheme of polynomials, so the surrogates can be expressed by a series of orthogonal polynomials w.r.t the distributions of the input variables. These polynomials can be extended as a basis of a polynomial space. In general, methods used to solve PCE problems are categorized as two types: intrusive and non-intrusive. The main idea behind the intrusive methods is the substitution of the input $\mathit{x}$ and output $f\left(\mathit{x}\right)$ with the truncated PCE and calculating the coefficients with the help of Galerkin projection [21]. However, the explicit formation of f is required to compose the Galerkin system, and a specific algorithm or program is needed to solve a particular problem. It is for these reasons that the intrusive models are not widely used; non-intrusive methods have been developed to avoid these limitations [21,22]. There are two main aspects of the non-intrusive methods: one is the choice of sampling strategies, for example Monte Carlo techniques; the other one is computational approaches. These two aspects are not independent of each other: for example, if $\mathit{x}\sim N(0,1)$, then the Gaussian quadrature method is introduced to solve the numerical integration and $\mathit{X}$ is the set of corresponding quadrature points. Another one of the more common methods in constructing surrogate models is the Gaussian process (GP), which is actually a Bayesian approach. Instead of attempting to identify a specific real model of the system, the GP method provides a posterior distribution over the model in order to make robust predictions about the system. As described in the highly influential works [23,24,25,26], the GP can be treated as a distribution over functions with properties controlled by a kernel. For the two prerequisites discussed in the previous paragraph, the GP generates a surrogate model that lies in a space spanned by kernels; meanwhile, Bayesian linear regression or classification methods are introduced to utilize the prior information.

Both the PCE and GP methods build surrogates, but there are some differences between them. The PCE method builds surrogates of a random variable $\mathit{y}$ as a function of another prior random variable $\mathit{x}$ rather than the distribution density function itself. The PCE surrogates are based on the orthogonal polynomial basis corresponding to the $p\left(\mathit{x}\right)$, so it is simple to obtain the mean and standard deviation of $\mathit{y}$. In contrast, the GP utilizes the covariance information so that it performs better in capturing the local features. Although both the PCE and GP approaches are feasible methods to compute the mean and standard deviation of $\mathit{y}$, the PCE performs more efficiently than the GP method.

As mentioned above, both the PCE and GP methods have their own trade-offs to consider when building surrogates, and there exists a connection to be explored. According to Paul Constantine’s work [27], ordinary kriging (i.e., GP in geostatistics) interpolation can be viewed as a transformed version of the least squares (LS) problem, and the PCE can be viewed as the least squares with selected basis and weights. However, the GP reverts to interpolation when the noise term is zero. When taking the noise term into consideration, the Gaussian process with the kernel (i.e., covariance matrix) $X^{T}X$ can be viewed as a ridge regression problem [28] with a regularization term. Furthermore, different numerical methods can affect the precision of the PCE method, as well. For example, Xiu [20] analyzed the aliasing error w.r.t the projection method and interpolation method. Thus, the inherent connection of the two models cannot be simply summarized as an LS solution, and how to output a model with high precision remains an interesting question.

There are connections between the PCE and GP methods that have been explored by R. Schobi, etc. They introduced a new meta-modeling method naming PC-kriging [29] (polynomial-chaos-based kriging) to solve the problems like rare event estimation [30], structural reliability analysis [31], quantile estimation [32], etc. In their papers, the PCE models can be viewed as a special form of GP where a Dirac function is introduced as the kernel. They also proposed the idea that the PCE models have better performance in capturing the global features and that the GP models approximate the local characteristics. We would like to describe the PC-kriging method as a GP model with a PCE-form trend function along with a noise term. The global features are dominated by the PCE trend, and local structures (residuals) are approximated by the ordinary GP process. The PC-kriging model thus introduces the coefficients as parameters to be optimized, and the solution can be derived by Bayesian linear regression with the basis consisting of the PCE polynomials. They also use the LARSalgorithms to calibrate the model and to select a sparse design. They construct a rigid framework to optimize the parameters, validate and calibrate the model and evaluate the model accuracy.

Unlike the PC-kriging, which takes the PCE as a trend, this paper focuses on the construction of the kernel in the GP to solve the regression problems, through which we can combine the two methods into a unified framework, unifying positive aspects from both and in so doing refining the surrogates. In other words, we wish to find the connection between the GP and the PCE by analyzing the attribution of their solutions, and we want to propose a new approach to achieve high-precision predictions. The main idea of this paper is described as follows. Firstly, the PCE surrogate is embedded in a Hilbert space whose bases are the orthonormal polynomials themselves, then a suitable inner product and a Mercer kernel [33] are defined to build a reproducing kernel Hilbert space (RKHS) [33]. Secondly, on the other hand, the kernel of the GP can be de-composited as the product of eigenfunctions, and we can define an inner product to generate a RKHS, as well. We have explicitly elaborated the two procedures respectively and proven that the two RKHS are isometrically isomorphic. Hence, a connection between these two approaches has been established via RKHS. Furthermore, we can obtain a solution of the PCE model by solving a GP model with the Mercer kernel w.r.t the PCE polynomial basis. We name this approach Gaussian processes on polynomial chaos basis (GPCB). In order to illustrate the capability of the GPCB method, we use the Kullback–Leibler divergence [34] to explicitly compare the PDFs of the posterior prediction of the GPCB and PCE method. Provided that the true function can be approximated by a finite number of PCE bases, it can be concluded that the GPCB can converge to the optimal subset of the RKHS wherein the true function lies.

The experimental design from the PCE model, i.e., the full tensor product of quadratures in each dimension, is used in the GPCB. We have overcome two concerns about the PCE and GP, respectively. Firstly, the PCE is based on a truncated polynomial basis, while the GPCB keeps all polynomials, which can be regarded as maintaining information in every feature. Secondly, the GP’s behavior depends on the experimental design; however, it often achieves the optimal result in local small datasets. The quadrature points derived from the PCE model are distributed evenly in the input space, and those points have high numerical precision w.r.t the polynomial basis; hence, they can work well with the GPCB. However, we must admit that the GPCB is still a GP approach, so when the dimension of input variables grows, the computational burden is on the table. In order to cope with the high dimensional problems, sampling strategies to lower the number of experimental designs are put on the table. The AK-MCSmethod [35] is a useful tool that adaptively selects new experimental designs; however, the experimental design tends to validate the selected surrogate model. We propose a new method that is model-free and that makes full use of the quadrature points. We randomly choose a sparse subset from the quadrature points to form a new experimental design while maintaining the accuracy. Several classical sampling strategies like MC, Halton and LHS are introduced to compare their capabilities. Our sample scheme has superior performance under the conditions in this paper. The GPCB is a novel method to build surrogate models, and it can be used for various physical problems such as reliability analysis and risk assessment.

This paper is divided into two parts. In Part 1, we discuss the mathematical rigor of the method: a brief summary of PCE and GP is presented in Section 2; the reproducing kernel Hilbert space (RKHS) is introduced to connect these two methods in Section 3; the GPCB method is proposed based on the discussion in Section 4; meanwhile, a theoretical Kullback–Leibler divergence between the GPCB and PCE method is demonstrated. In Part 2, an explicit Mehler kernel is presented with the Hermite polynomial basis in the last part of Section 4; several tests of the GPCB with some benchmark functions are presented in Section 5, along with the random constructive sampling method for high dimensional problems.

2. Brief Review of PCE and GP

Firstly, we want to have a clear idea of how the PCE and GP work under the circumstances that, e.g., only samples of input X and output Y are obtained. Different assumptions are made to cope with PCE and GP, respectively, and the processing procedures are presented in the following subsections.

2.1. Polynomial Chaos Expansion

Just as discussed in Section 2.2, the output is assumed to be represented by a model $y=f\left(\mathit{x}\right)+ϵ,$ where $f\left(\mathit{x}\right):\mathit{x}\in \Omega \to \mathbb{R}$ is the real function underneath and $ϵ\sim \mathcal{N}(0,{\sigma }_{ϵ}^2)$. Here, we define $\mathit{x}={(x_1,x_2,\dots ,x_d)}^T$ as a d-dimensional vector of independent random variables in a bounded domain $\Omega \subset {\mathbb{R}}^d$. Suppose $\{x_i,i=1,\dots ,d\}$ are independent and identically distributed; the joint PDF has the form $p\left(\mathit{x}\right)={\prod }_{i=1}^{d}p\left(x_i\right)$. In the context of PCE, we aim to seek a surrogate of the model $f\left(\mathit{x}\right)$ as an expansion of a series of orthonormal polynomials ${\mathbf{\varphi }}_{\alpha }\left(\mathit{x}\right)$:

$$\begin{array}{c}\hfill f\left(\mathit{x}\right)=\sum _{\alpha \in {\mathbb{N}}^d}{\beta }_{\alpha }{\mathbf{\varphi }}_{\alpha }\left(\mathit{x}\right),\alpha =\{{\alpha }_1,\dots ,{\alpha }_d\}\end{array}$$

(1)

where $\alpha$ is the multi-index, ${\mathbf{\varphi }}_{\alpha }\left(\mathit{x}\right)={\prod }_{i=1}^d{\varphi }_{{\alpha }_i}^{\left(i\right)}\left(x_i\right)$, ${\int }_{{\Omega }_i}{\varphi }_m^{\left(i\right)}\left(x_i\right){\varphi }_n^{\left(i\right)}\left(x_i\right)p_i\left(x_i\right)d{x}_i={\delta }_{mn}$ with ${\Omega }_i$ the marginal domain of $\Omega$ and ${\delta }_{mn}$ the Kronecker delta. Xiu et al. [20] have summarized various correspondences between the distribution and polynomial basis to form generalized polynomial chaos.

It is proven that the original model $f\left(\mathit{x}\right)$ can be approximated to any degree of accuracy in a strong sense [20], e.g., mean-square norm $\parallel f\left(\mathit{x}\right)-{\sum }_{\alpha \in {\mathbb{N}}^D}{\beta }_{\alpha }{\mathbf{\varphi }}_{\alpha }\left(\mathit{x}\right)\parallel$ in an $L^2$ norm defined on $\Omega$, although f is not necessarily the span of orthonormal polynomial bases. Since we are unable to calculate an infinite series, the truncation scheme corresponding to multi-index $\alpha$ is introduced such that we can rearrange the polynomials. For simplicity, we can rewrite Equation (1) in the following form:

$$\begin{array}{c}\hfill f\left(\mathit{x}\right)\approxeq \sum _{l=0}^M{\beta }_l{\mathbf{\varphi }}_l\left(\mathit{x}\right)\triangleq f_P,\end{array}$$

(2)

We can simply solve the above system via the ordinary least squares method or the non-intrusive method. Specifically, we focus on the non-intrusive projection method, whereby we can directly obtain the coefficients by taking the expectation value of Equation (2) multiplied by ${\mathbf{\varphi }}_l\left(\mathit{x}\right)$:

$$\begin{array}{c}\hfill {\beta }_l=\int f\left(\mathit{x}\right){\mathbf{\varphi }}_l\left(\mathit{x}\right)p\left(\mathit{x}\right)d\mathit{x}\approx \sum _{i=1}^N{\omega }_{i}f\left(X_i\right){\mathbf{\varphi }}_l\left(X_i\right),l=0,\dots ,M\end{array}$$

(3)

where the second equation is derived by the numerical integration techniques, such as the Gaussian quadrature rule, and $\{X_i,i=1,\dots ,N\}$ and $\{{\omega }_i,i=1,\dots ,N\}$ are the corresponding nodes and weights. The integration is exact when $f\left(\mathit{x}\right)$ is of polynomial complexity. Together with Equations (2) and (3)$,f_P\left(\mathit{x}\right)$ has the form:

$$\begin{array}{c}\hfill f_P\left(\mathit{x}\right)\approx \sum _{l=0}^M\left[\sum _{i=1}^N{\omega }_{i}f\left(X_i\right){\mathbf{\varphi }}_l\left(X_i\right)\right]{\mathbf{\varphi }}_l\left(\mathit{x}\right)\triangleq \sum _{l=0}^M{\beta }_l{\mathbf{\varphi }}_l\left(\mathit{x}\right).\end{array}$$

(4)

$\{f\left(X_i\right),i=1,\dots ,N\}$ remain unknown to us, and usually, they are substituted by $\{Y_i,i=1,\dots ,N\}$. Note $Y_i=f\left(X_i\right)+{ϵ}_i$, so such a substitution will introduce noise into the surrogate; hence, the approximation error is neglected as a source of uncertainty.

2.2. Gaussian Process Regression

The analysis of the Gaussian process regression model [26] is reviewed in this section. A Gaussian prior is placed over function $f\left(\mathit{x}\right)$, i.e., $f\left(\mathit{x}\right)\sim \mathcal{GP}\left(m\right(\mathit{x}),k(\mathit{x},\mathit{x}\left)\right)$, where m is the mean function and k is the kernel function, which is positive semi-definite bounded. More specifically, let $\mathit{X}=\{X_i,i=1,\dots ,N\}\in {\Omega }^N$ be the input data, and let $Y=\{Y_i,i=1,\dots ,N\}\in {\mathbb{R}}^N$ be the output data, then we have $Y=f\left(\mathit{X}\right)+\mathit{ϵ}$ with $f\left(\mathit{X}\right)\sim \mathcal{N}\left(m\right(\mathit{X}),k(\mathit{X},\mathit{X}\left)\right)$ and $\mathit{ϵ}\sim \mathcal{N}(0,{\sigma }_{ϵ}^{2}I)$. Bear in mind that the mathematical expression of $f\left(\mathit{x}\right)$ is implicit, so $f\left(\mathit{x}\right)$ is approximated to achieve the best guess prediction $f_G$ in the statistical sense. With the help of Bayes’ theorem, prediction and corresponding variance at a new point $\mathit{x}$ can be obtained by the following equations [36]:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_G\left(f\left(\mathit{x}\right)\right|Y,\mathit{X},\mathit{x},\mathbf{\theta })& =\mathcal{N}(f_G\left(\mathit{x}\right),cov\left(f_G\left(\mathit{x}\right)\right)),\hfill \\ \hfill f_G\left(\mathit{x}\right)\triangleq \mathbb{E}\left[f\right(\mathit{x}\left)\right|Y,\mathit{X},\mathit{x},\mathbf{\theta }]& =K_{\mathit{x}}^T{[K+{\sigma }_{ϵ}^{2}I]}^{-1}Y,\hfill \\ \hfill cov\left(f_G\left(\mathit{x}\right)\right)& =K_{\mathbf{xx}}-K_{\mathit{x}}^T{[K+{\sigma }_{ϵ}^{2}I]}^{-1}{K}_{\mathit{x}}.\hfill \end{array}\end{array}$$

(5)

where $K=k(\mathit{X},\mathit{X})\in {\mathbb{R}}^{N\times N}$ as the covariance matrix with $K_{ij}=k(X_i,X_j)$ and $K_{\mathit{x}}=k(\mathit{X},\mathit{x})\in {\mathbb{R}}^{N\times 1}$, $K_{\mathit{xx}}=k(\mathit{x},\mathit{x})\in \mathbb{R}$ are defined similarly. Note that Equation (5) shows that the mean value of the posterior distribution can be expressed as a linear combination of N kernel functions as follows:

$$\begin{array}{c}\hfill f_G\left(\mathit{x}\right)=\sum _{i=1}^N{\alpha }_{i}k(X_i,\mathit{x}),\alpha ={(K+{\sigma }_{ϵ}^{2}I)}^{-1}Y\end{array}$$

(6)

3. Links between the PCE and GP

The basic concepts of PCE and GP are discussed in Section 2. GP generates a surrogate based on Bayes’ theorem and the Gaussian hypothesis; however, it is controlled by the kernel function and the experimental design and usually does not utilize prior distribution information; PCE substitutes the model with orthonormal polynomials, which is more computational efficient, but performs badly when facing noisy or big data. This section aims to build a connection between PCE and GP; hence, they can be studied in the same structure and be combined to improve the performance of the surrogates. The reproducing kernel Hilbert space will be of great help to build such a bridge, and we are going to present it as follows.

3.1. Generate an RKHS from a Mercer Kernel Constructed by the PCE Basis

We have obtained a complete orthonormal basis $\left\{{\mathbf{\varphi }}_l\right\}$ of Hilbert space $\mathcal{H}:=span\left\{{\mathbf{\varphi }}_l\left(\mathit{x}\right)\right\}$ with inner product $<f\left(\mathit{x}\right),g\left(\mathit{x}\right)>=\int f\left(\mathit{x}\right)g\left(\mathit{x}\right)p\left(\mathit{x}\right)d\mathit{x}$ in Section 2.1. We can see that the PCE method generates surrogates, which are actually a linear combination of $\left\{{\mathbf{\varphi }}_l\right\}$, so there exists a unique expansion $f={\sum }_l{f}_l{\mathbf{\varphi }}_l\in \mathcal{H}$. According to Mercer’s theorem [33], we aim to define a kernel having the following form:

$$\begin{array}{c}\hfill k(\mathit{x},{\mathit{x}}^{\prime })=\sum _l{\lambda }_l{\mathbf{\varphi }}_l\left(\mathit{x}\right){\mathbf{\varphi }}_l\left({\mathit{x}}^{\prime }\right)s.t.k(\mathit{x},\mathit{x})<\infty for\forall \mathit{x}\in \Omega .\end{array}$$

(7)

If we have positive weights ${\lambda }_l$ that satisfy ${\sum }_l{\lambda }_l{\mathbf{\varphi }}_l^2\left(\mathit{x}\right)<\infty$, then for any $\mathit{x}\in \Omega$, together with the Cauchy–Schwarz inequality, we have:

$$\begin{array}{c}\hfill \left|f\left(\mathit{x}\right)\right|\le \sqrt{\sum _l\frac{<f\left(\mathit{x}\right),{\mathbf{\varphi }}_l\left(\mathit{x}\right){>}^2}{{\lambda }_l}}\sqrt{\sum _l{\lambda }_l{\mathbf{\varphi }}_l^2\left(\mathit{x}\right)}.\end{array}$$

(8)

$\left|f\right(\mathit{x}\left)\right|$ is point-wise bounded because $f\left(\mathit{x}\right)\in \mathcal{H}$ for any $\mathit{x}\in \Omega$. By checking the right side of the above inequality, the second term is ensured in advance, then $f\left(\mathit{x}\right)$ lies in a subspace of $\mathcal{H}$ such that:

$$\begin{array}{c}\hfill {\mathcal{H}}_P=\left\{f\in \mathcal{H}|<f,f{>}_{{\mathcal{H}}_P}=\sum _l\frac{<f\left(\mathit{x}\right),{\mathbf{\varphi }}_l\left(\mathit{x}\right){>}^2}{{\lambda }_l}<\infty ,\sum _l{\lambda }_l{\mathbf{\varphi }}_l^2\left(\mathit{x}\right)<\infty \right\}.\end{array}$$

(9)

Proposition 1.

${\mathcal{H}}_P$ defined in Equation (9) is an RKHS with Mercer kernel defined in Equation (7).

3.2. Generate an RKHS from the Reproducing Kernel Map Construction

We aim to compose a space of functions in which all the GP surrogates are embedded. Given Equation (6) and an arbitrary experimental design $\mathit{X}$, define a space of functions as follows:

$$\begin{array}{c}\hfill {\mathcal{H}}_G^{\prime }=\left\{f\left(\mathit{x}\right)=\sum _{i=1}^N{f}_{i}k(\mathit{x},X_i)|N\in \mathbb{N},\mathit{X}\in {\Omega }^N,\mathit{x}\in \Omega ,f_i\in \mathbb{R},\sum _{i=1}^N\sum _{j=1}^N{f}_i{f}_{j}k(X_i,X_j)<+\infty \right\}.\end{array}$$

(10)

Proposition 2.

${\mathcal{H}}_G^{\prime }$ is a pre-Hilbert space with the inner product $<·,·{>}_{{\mathcal{H}}_G^{\prime }}$

Now that ${\mathcal{H}}_G^{\prime }$ is a pre-Hilbert space and given the norm ${\parallel f\left(\mathit{x}\right)\parallel }_{{\mathcal{H}}_G^{\prime }}=\sqrt{<f\left(\mathit{x}\right),f\left(\mathit{x}\right){>}_{{\mathcal{H}}_G^{\prime }}}$, we can define a closure of ${\mathcal{H}}_G^{\prime }$ as ${\mathcal{H}}_G$ derived by the classical Hilbert space theory. This is an abstract space where the norm of ${\mathcal{H}}_G^{\prime }$ extends to the closure ${\mathcal{H}}_G$. Thus, we have a Hilbert space ${\mathcal{H}}_G$.

Proposition 3.

${\mathcal{H}}_G$ defined above is the unique RKHS of the kernel $k(·,·)$.

3.3. Reproducing Kernel Hilbert Spaces as a Linkage

${\mathcal{H}}_P$ and ${\mathcal{H}}_G$ are RKHS with the Mercer kernel and GP kernel, respectively. We are going to investigate the relationship between the two RKHS, by which we can discuss the two approaches in a unified structure. Let $\mathit{X}$ be a sample set and GP kernel $k(·,·)$ be a real positive semi-definite kernel, then according to Mercer’s theorem, $k(X_i,X_j)$ has an eigenfunction expansion:

$$\begin{array}{c}\hfill k(X_i,X_j)=\sum _l{\lambda }_l{\mathbf{\varphi }}_l\left(X_i\right){\mathbf{\varphi }}_l\left(X_j\right),\end{array}$$

(11)

where the eigenfunctions $\left\{{\mathbf{\varphi }}_l\right\}$ are orthonormal, i.e., $<{\mathbf{\varphi }}_l,{\mathbf{\varphi }}_{l^{\prime }}>=\int {\mathbf{\varphi }}_l\left(\mathit{x}\right){\mathbf{\varphi }}_{l^{\prime }}\left(\mathit{x}\right)p\left(\mathit{x}\right)d\mathit{x}={\delta }_{l{l}^{\prime }}$ and $\{{\lambda }_l,{\mathbf{\varphi }}_l\}$ satisfies ${\sum }_l{\lambda }_l{\mathbf{\varphi }}_l^2\left(\mathit{x}\right)<\infty$. Let $f_{\mathit{X}}\left(\mathit{x}\right)\in {\mathcal{H}}_G$ with experimental design $\mathit{X}$, then we can rewrite it according to Equations (10) and (11):

$$\begin{array}{c}\hfill f_{\mathit{X}}\left(\mathit{x}\right)=\sum _{i=1}^N{f}_i\sum _l{\lambda }_l{\mathbf{\varphi }}_l\left(X_i\right){\mathbf{\varphi }}_l\left(\mathit{x}\right)=\sum _l\left[{\lambda }_l\sum _{i=1}^N{f}_i{\mathbf{\varphi }}_l\left(X_i\right)\right]{\mathbf{\varphi }}_l\left(\mathit{x}\right)\triangleq \sum _l{c}_l\left(\mathit{X}\right){\mathbf{\varphi }}_l\left(\mathit{x}\right).\end{array}$$

(12)

where $c_l\left(\mathit{X}\right)$ is identically determined by l and $\mathit{X}$, and it has a similar form as a function lies in ${\mathcal{H}}_P$. Actually, given $f_{\mathit{X}}\left(\mathit{x}\right),g_{{\mathit{X}}^{\prime }}\left(\mathit{x}\right)\in {\mathcal{H}}_G$, we have:

$$\begin{array}{cc}\hfill <f_{\mathit{X}}\left(\mathit{x}\right),g_{{\mathit{X}}^{\prime }}\left(\mathit{x}\right){>}_{{\mathcal{H}}_G}& =\sum _{i=1}^N\sum _{j=1}^{N^{\prime }}{f}_i{g}_j\left[\sum _l{\lambda }_l{\mathbf{\varphi }}_l\left(X_i\right){\lambda }_l{\mathbf{\varphi }}_l\left(X_j^{\prime }\right)/{\lambda }_l\right]\hfill \\ & =\sum _i\left\{\left[{\lambda }_l\sum _{i=1}^N{f}_i{\mathbf{\varphi }}_l\left(X_i\right)\right]\left[{\lambda }_l\sum _{j=1}^{N^{\prime }}{g}_j{\mathbf{\varphi }}_l\left(X_j^{\prime }\right)\right]\right\}/{\lambda }_l\hfill \\ & =\sum _l{c}_l\left(\mathit{X}\right)c_l\left({\mathit{X}}^{\prime }\right)/{\lambda }_l\hfill \\ & =\sum _l<f_{\mathit{X}}\left(\mathit{x}\right),{\mathbf{\varphi }}_l><g_{{\mathit{X}}^{\prime }}\left(\mathit{x}\right),{\mathbf{\varphi }}_l>/{\lambda }_l\hfill \\ & =<f_{\mathit{X}}\left(\mathit{x}\right),g_{{\mathit{X}}^{\prime }}\left(\mathit{x}\right){>}_{{\mathcal{H}}_P}.\hfill \end{array}$$

(13)

The above equation gives us the information that their inner product stands in ${\mathcal{H}}_P$, as well, so we can conclude that $f_{\mathit{X}}$ lies in ${\mathcal{H}}_P$. It also shows us that the two inner products are equivalent. Next, we are going to propose a rigid theorem to prove that the two spaces are isometrically isomorphic.

Theorem 1.

The reproducing kernel Hilbert space ${\mathcal{H}}_G$ of a given kernel k is isometrically isomorphic to the space ${\mathcal{H}}_P$.

According to the proof of Theorem 1 in Appendix A.4, it is reasonable to introduce a weighted $l^2$ space $l_{1/\lambda }^2$ because it is difficult to find a direct linear map between ${\mathcal{H}}_G$ (where $c_l$ varies according to $\mathit{X}$ and k) and ${\mathcal{H}}_P$ (where $c_l$ varies according to the distribution of $\mathit{x}$). Figure 1 shows two flowcharts used to describe different processes in generating the RKHS.

Furthermore, the GP prediction is a combination of the kernel functions, which consist of infinite eigenfunctions, while the PCE prediction is always a combination of finite polynomial bases. The Kullback–Leibler divergence (KL divergence) is a useful criterion to indicate the performance of different surrogate models. We are going to present the comparison of the GPCB and PCE methods with the help of KL divergence in the next section.

4. Gaussian Process on Polynomial Chaos Basis

${\mathcal{H}}_G$ and ${\mathcal{H}}_P$ are isomorphic as discussed in previous Section 3, so it is natural to come up with the idea that GP can be conducted with $k(·,·)$ as the Mercer Kernel generated by polynomial basis in the PCE, and the new model is called Gaussian process on polynomial chaos basis (GPCB). In fact, the GPCB generates a PCE-like model, but with a different philosophy. Note that the posterior distribution of the predictions regarding experimental design $\{\mathit{X},Y\}$ can be calculated analytically, so we are able to compute the KL divergence as well, which are presented as follows.

4.1. Comparison of the PCE and GPCB with the Kullback–Leibler Divergence

The true distribution of the system is always implicit in practice. Without loss of generality, the underlying true system is assumed to be $f_{\overline{P}}\left(\mathit{x}\right)={\sum }_{l=0}^{\overline{M}}{\overline{\beta }}_l{\mathbf{\varphi }}_l\left(\mathit{x}\right)$ if $f_{\overline{P}}\left(\mathit{x}\right)\in C^0\left(\overline{\Omega }\right)$ such that it can be approximated by the polynomials to any degree of accuracy [20].

Firstly, we presume that $\overline{M}\le M$, i.e., $\mathbf{\beta }$ in Equation (4) is an unbiased estimator of $\overline{\mathbf{\beta }}$. Hence, the PCE approximation can be considered as a precise approximate of the true function. We compare the performance of the GPCB and the PCE method by comparing their difference in the posterior distribution of the prediction. It is known that given experimental design $\{\mathit{X},Y\}$ and kernel function $k(\mathit{x},{\mathit{x}}^{\prime })={\sum }_{l=0}^{\infty }{\lambda }_l{\mathbf{\varphi }}_l\left(\mathit{x}\right){\mathbf{\varphi }}_l\left({\mathit{x}}^{\prime }\right)$, the distribution of the prediction of the GPCB reads:

$$\begin{array}{c}\hfill p_G\left(f_G\left(\mathit{x}\right)\right)=\mathcal{N}(K_{\mathit{x}}^T{K}_Y^{-1}Y,K_{\mathit{xx}}-K_{\mathit{x}}^T{K}_Y^{-1}{K}_{\mathit{x}})\triangleq \mathcal{N}({\mu }_1,{\Sigma }_1),\end{array}$$

(14)

where conditions $\mathit{X},Y,\mathit{x}$ are dropped in $p_G\left(f_G\left(\mathit{x}\right)\right|\mathit{X},Y,\mathit{x})$ for simplicity and $K_Y=K+{\sigma }^{2}I$. Similarly, the prediction of the PCE with the projection method is $f_P\left(\mathit{x}\right)=\mathbf{\varphi }{\Phi }^{T}WY$, which is derived from the estimation $\mathbf{\beta }={\Phi }^{T}WY$. Here, $\mathbf{\varphi }=\mathbf{\varphi }\left(\mathit{x}\right)\in {\mathbb{R}}^{1\times (P+1)}$, $\Phi =\mathbf{\varphi }\left(\mathit{X}\right)\in {\mathbb{R}}^{N\times (P+1)}$, $W=diag\{{\omega }_1,\dots ,{\omega }_N\}$ is a diagonal matrix. The corresponding prediction variance is $cov\left(f_P\left(\mathit{x}\right)\right)=\mathbf{\varphi }cov\left(\mathbf{\beta }\right){\mathbf{\varphi }}^T={\sigma }^2\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T$. Dropping the conditions in $p_P\left(f_P\left(\mathit{x}\right)\right|\mathit{X},Y,\mathit{x})$ as well, the previous results indicate:

$$\begin{array}{c}\hfill p_P\left(f_P\left(\mathit{x}\right)\right)=\mathcal{N}(\mathbf{\varphi }{\Phi }^{T}WY,{\sigma }^2\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T)\triangleq \mathcal{N}({\mu }_2,{\Sigma }_2).\end{array}$$

(15)

We can evaluate the discrepancy between $p_G\left(f_G\left(\mathit{x}\right)\right)$ and $p_P\left(f_P\left(\mathit{x}\right)\right)$, hence comparing their performance. The KL divergence can be calculated analytically:

$$\begin{array}{c}\hfill D_{KL}\left(p_P,p_G\right)=\frac{1}{2}\left(-1+\frac{{\Sigma }_2}{{\Sigma }_1}+log\frac{{\Sigma }_1}{{\Sigma }_2}+\frac{{\left({\mu }_1-{\mu }_2\right)}^2}{{\Sigma }_1}\right)\triangleq \frac{1}{2}\left(-1+b-logb+\frac{a^2}{{\Sigma }_2}b\right),\end{array}$$

(16)

where $a,b$ are simplified notations for the corresponding parts in Equation (16). In fact, b is the point-wise ratio between the posterior variances of the predictions, and a represents the difference between the posterior mean of the prediction. We discuss the properties of $D_{KL}$ starting from a special case to the general conditions hereafter.

Let k be a truncated kernel with the M basis of PCE, i.e., $k(\mathit{x},{\mathit{x}}^{\prime })={\sum }_{l=0}^M{\lambda }_l{\mathbf{\varphi }}_l\left(\mathit{x}\right){\mathbf{\varphi }}_l\left({\mathit{x}}^{\prime }\right)$. Actually, it can be seen as the assignment of $\{{\lambda }_l,l>M\}$ with the value of zero, which can be achieved by optimizing the value ${\lambda }_l$ with a specific procedure. b can be simplified as:

$$\begin{array}{cc}\hfill b& =\frac{{\sigma }^2\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T}{\mathbf{\varphi }\left(\Lambda -\Lambda {\Phi }^T{K}_Y^{-1}\Phi \Lambda \right){\mathbf{\varphi }}^T}=\frac{\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T}{\mathbf{\varphi }{\Phi }^{T}W{K}_Y^{-1}KW\Phi {\mathbf{\varphi }}^T}\hfill \\ & =\frac{\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T}{\mathbf{\varphi }{\Phi }^{T}WU{(S+{\sigma }^{2}I)}^{-1}{U}^{T}US{U}^{T}W\Phi {\mathbf{\varphi }}^T}\in \left[\frac{s_{max}+{\sigma }_{ϵ}^2}{s_{max}},\frac{s_{min}+{\sigma }_{ϵ}^2}{s_{min}}\right],\hfill \end{array}$$

(17)

where $\Lambda =diag\left\{{\lambda }_l\right\}$ is a diagonal matrix and $K=US{U}^T$ is the eigenvalue decomposition. Let $s_{max}$ be the maximum eigenvalue and $s_{min}$ be the minimal one; the above interval holds because K is a positive definite matrix. Note that b is an invariant with fixed $\mathit{x}$. On the other hand, the distribution of a is as follows:

$$\begin{array}{cc}\hfill a& =\varphi \left({\Phi }^{T}WY-\Lambda {\Phi }^T{K}_Y^{-1}Y\right)={\sigma }_{ϵ}^2\mathbf{\varphi }{\Phi }^{T}W{K}_Y^{-1}Y\hfill \\ & \sim \mathcal{N}\left({\sigma }_{ϵ}^2\mathbf{\varphi }{\Phi }^{T}W{K}_Y^{-1}\widehat{Y},{\sigma }_{ϵ}^6\mathbf{\varphi }{\Phi }^{T}W{K}_Y^{-1}{K}_Y^{-1}W\Phi {\mathbf{\varphi }}^T\right).\hfill \end{array}$$

(18)

Here, $\widehat{Y}$ denotes the mean value of the observations, i.e., the true response. It is necessary to state that the randomness of a is brought by the random variable $\mathit{ϵ}$ in observation Y. In fact, we have the expectation of $D_{KL}(p_P,p_G)$ as follows:

$$\begin{array}{cc}\hfill {\mathbb{E}}_{ϵ}\left[D_{KL}(p_P,p_G)\right]& =\frac{1}{2}\left(-1+b-logb+\frac{{\mathbb{E}}_{ϵ}{a}^2}{{\Sigma }_2}b\right)\hfill \\ & =\frac{1}{2}\left(-1+b-logb+(var\left(c\right)+{\left({\mathbb{E}}_{ϵ}a\right)}^2)b/{\Sigma }_2\right)\hfill \\ & \le \frac{1}{2}\left(-1+b-logb+\frac{{\sigma }_{ϵ}^2+{\parallel \widehat{Y}\parallel }^2}{{\sigma }_{ϵ}^2}{\left(\frac{{\sigma }_{ϵ}^2}{s_{min}+{\sigma }_{ϵ}^2}\right)}^{2}b\right).\hfill \end{array}$$

(19)

Presume that the observation Y is normalized, as well as $\widehat{Y}$; hence, ${\parallel \widehat{Y}\parallel }^2$ can be estimated as $\mathcal{O}\left(1\right)$. The main difference is affected by the kernel $k(·,·)$ (or $\Lambda$) and the term ${\sigma }_{ϵ}^2$. More specifically, The GPCB can achieve a smaller variance than the PCE method in a point-wise manner because $b>1$, and the differences between the predictions of the two methods is of the order of ${\sigma }_{ϵ}^2$. Furthermore, if ${\sigma }_{ϵ}^2$ is sufficient small, i.e., ${\sigma }_{ϵ}^2\ll s_{min}$, we have $b\to 1$, thus ${\mathbb{E}}_{ϵ}\left[D_{KL}(p_P,p_G)\right]\to 0$. If ${\sigma }_{ϵ}^2=0$, i.e., we investigate the noise-free models, then $b=1$ and $a=0$, which enforces ${\Sigma }_1={\Sigma }_2$ and ${\mu }_1={\mu }_2$, respectively. This means that $p_G$ and $p_P$ are identical distributions, i.e., $D_{KL}\left(p_P,p_G\right)=0$.

We can conclude that the expected value of $D_{KL}(p_P,p_G)$ is bounded by a certain constant, which mainly depends on the ${\sigma }_{ϵ}^2$ and $\Lambda$. In other words, since ${\sigma }_{ϵ}^2$ is given in the prior, and the $\Lambda$ are optimized; hence, the $D_{KL}(p_P,p_G)$ is constrained, which means that the GPCB is as stable as the PCE method. Nonetheless, if $f_P$ has reached a desired prediction precision, then $f_G$ with the kernel constructed with the same basis can have a desirable precision, as well as smaller variance.

Secondly, we consider that $\overline{M}>M$, where the $\mathbf{\beta }$ of the PCE method is not an unbiased estimate of $\overline{\mathbf{\beta }}$ any longer. Let $p_{\overline{P}}$ denote the PCE approximation with the $\overline{M}$ basis; under the circumstance, $k(\mathit{x},{\mathit{x}}^{\prime })={\sum }_{k=0}^{\overline{M}}{\lambda }_l{\mathbf{\varphi }}_l\left(\mathit{x}\right){\mathbf{\varphi }}_l\left({\mathit{x}}^{\prime }\right)$ is achieved by tuning the value of $\Lambda$ via a certain learning method; hence, $D_{KL}(p_{\overline{P}},p_G)$ is also bounded by a constant according to Equation (19), i.e., the GPCB can converge to the precise PCE prediction $p_{\overline{P}}$, as well. However, the KL divergence $D_{KL}(p_{\overline{P}},p_P)$ is given as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_{KL}\left(p_{\overline{P}},p_P\right)& =\frac{1}{2}\left(-1+\frac{\overline{\Sigma }}{{\Sigma }_2}-log\frac{\overline{\Sigma }}{{\Sigma }_2}+\frac{{\left(\overline{\mu }-{\mu }_2\right)}^2}{{\Sigma }_1}\right)\triangleq \frac{1}{2}\left(-1+\overline{b}-log\overline{b}+\frac{{\overline{a}}^2}{{\Sigma }_1}\right),\hfill \\ \hfill \overline{b}& =\frac{\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T+{\mathbf{\varphi }}_r{\Phi }_r^T{W}^2{\Phi }_r{\mathbf{\varphi }}_r^T+\mathbf{\varphi }A{\mathbf{\varphi }}_r^T}{\mathbf{\varphi }{\Phi }^T{W}^2\Phi {\mathbf{\varphi }}^T}\ge 1,\hfill \\ \hfill \overline{a}& ={\mathbf{\varphi }}_r{\Phi }_r^{T}WY\sim \mathcal{N}\left({\mathbf{\varphi }}_r{\Phi }_r^{T}W\widehat{Y},{\sigma }_{ϵ}^2{\mathbf{\varphi }}_r{\Phi }_r^T{W}^2{\Phi }_r{\mathbf{\varphi }}_r^T\right).\hfill \end{array}\end{array}$$

(20)

We denote ${\mathbf{\varphi }}_r$ as the basis that belongs to the model $f_{\overline{P}}$, but $f_P$. It is shown that the biased PCE $f_P$ has smaller variance, however with a bias whose mean value depends on ${\mathbf{\varphi }}_r$. We notice that ${\mathbf{\varphi }}_r$ represents the high-order polynomials; hence, the bias can be considerably large for general cases, and so is the $D_{KL}\left(p_{\overline{P}},p_P\right)$. We can conclude that even though the biased PCE has smaller variance, the relatively large bias can lead to a false prediction.

In fact, if the underlying system is smooth enough to be modeled by a polynomial approximation, then we can adaptively increase the number of polynomial bases (and the experimental designs if necessary) to reach a precise approximation. However, on the other hand, we can directly use the GPCB method, which is a one-step Bayesian approximation, that converges to the hypothetical true system $f_{\overline{P}}$. Roughly speaking, the GPCB finds $\overline{M}$ automatically by tuning the parameters $\Lambda$ instead of adaptively changing the value of P in the PCE method. It indeed provides more convenience for computation. The key problem is the evaluation of $\Lambda$. Specifically, we introduce the Mehler kernel [37], which is an analytic expression of the Mercer kernel constructed by Hermite polynomials, and discuss the learning procedure of ${\Lambda }_l$.

4.2. Construction of the Kernel with Hermite Polynomials

Recall that we have $\left\{{\mathbf{\varphi }}_{\alpha }\right\}$ in PCE as an orthonormal basis, then we regard them as eigenfunctions of a kernel $k(·,·)$. Since $p\left(\mathit{x}\right)$ follows the standard Gaussian distribution, then ${\varphi }_l^{\left(i\right)}\left(x_i\right)=H{e}_l\left(x_i\right)/\sqrt{l!}$, where $x_i$ is the i-th variable of $\mathit{x}$ and $H{e}_l\left(x_i\right)$ is the Hermite polynomial of degree l:

$$\begin{array}{c}\hfill H{e}_l\left(x_i\right)={(-1)}^l{e}^{\frac{x_i^2}{2}}\frac{d^l}{d{x}_i^l}{e}^{-\frac{x_i^2}{2}}.\end{array}$$

(21)

Here, we denote the real l-th multi-index of the l-th polynomial in Equation (2) as ${\mathbf{\alpha }}^{\left(l\right)}=({\alpha }_1^{\left(l\right)},\dots ,{\alpha }_d^{\left(l\right)})$, $|{\mathbf{\alpha }}^{\left(l\right)}|={\alpha }_1^{\left(l\right)}+\dots +{\alpha }_d^{\left(l\right)}$ and ${\mathbf{\alpha }}^{\left(l\right)}!={\alpha }_1^{\left(l\right)}!\dots {\alpha }_d^{\left(l\right)}!$, then according to the previous analysis, we can get ${\mathbf{\varphi }}_{{\mathbf{\alpha }}^{\left(l\right)}}\left(\mathit{x}\right)={\mathit{He}}_{{\mathbf{\alpha }}^{\left(l\right)}}\left(\mathit{x}\right)/\sqrt{{\mathbf{\alpha }}^{\left(l\right)}!}$. We have Mehler kernel $Me(\mathit{X},{\mathit{X}}^{\prime })$ [37] with orthonormal Hermite polynomials as eigenfunctions:

$$\begin{array}{c}\hfill Me(\mathit{x},{\mathit{x}}^{\prime })=exp\left\{-\frac{D_{\mathit{x}}{D}_{\mathit{x}}^T+D_{{\mathit{x}}^{\prime }}{D}_{{\mathit{x}}^{\prime }}^T}{2}\right\}\left[exp\sum _{i=1}^d{\rho }_i{x}_i{x}_i^{\prime }\right]=\sum _{l=0}^{\infty }{\mathbf{\rho }}^{{\mathbf{\alpha }}^{\left(l\right)}}{\mathbf{\varphi }}_{{\mathbf{\alpha }}^{\left(l\right)}}\left(\mathit{x}\right){\mathbf{\varphi }}_{{\mathbf{\alpha }}^{\left(l\right)}}\left({\mathit{x}}^{\prime }\right),\end{array}$$

(22)

where the eigenvalue is ${\lambda }_{{\mathbf{\alpha }}^{\left(l\right)}}={\mathbf{\rho }}^{{\mathbf{\alpha }}^{\left(l\right)}}={\prod }_{i=1}^d{\rho }_i^{{\alpha }_i^{\left(l\right)}}$ for parameter $\mathbf{\rho }$, $D_x$ is the symbol representing the row gradient operator, i.e., $D_{\mathit{x}}=(\partial /\partial x_1,\dots ,\partial /\partial x_d)$, and $D_{{\mathit{x}}^{\prime }}$ is defined similarly. Specifically, in the one-dimensional case:

$$\begin{array}{c}\hfill Me(x,x^{\prime })=\frac{1}{\sqrt{1-{\rho }^2}}exp\left(-\frac{{\rho }^2(x^2+{x\prime }^2)-2\rho x{x}^{\prime }}{2(1-{\rho }^2)}\right)=\sum _{l=0}^{\infty }{\rho }^l{\varphi }_l\left(x\right){\varphi }_l\left(x^{\prime }\right),\end{array}$$

(23)

where the eigenvalue ${\lambda }_l={\rho }^l>0$. The truncated kernel $M{e}_M(x,x^{\prime })={\sum }_{l=0}^M{\rho }^l{\varphi }_l\left(x\right){\varphi }_l\left(x^{\prime }\right)$, and its attribution can be investigated by varying $\rho$ and M. Figure 2a illustrates the truncated kernel $M{e}_M(x,x^{\prime })$, which shows that $M{e}_M(x,x^{\prime })$ tends to converge to $Me(x,x^{\prime })$ as M grows. Figure 2b shows the values of $Me(x,-0.8)$ with different $\rho$. It presents to us that the influence of eigenvalue ${\lambda }_l$ is greater on the Mehler kernel.

4.3. Learning the Hyper-Parameter $\rho$ of $Me(x,x^{\prime })$

It is clear that ${\lambda }_l={\rho }^l$ has a great impact on the kernel values, hence affecting the convergence of the $f_G\left(\mathit{x}\right)$. We start with a simple example, where $f\left(x\right)=5+x+exp\left(x\right),x\sim \mathcal{N}(0,2^2)$ is the true underlying function, and the noise term is ignored in the observation. Considering the Taylor expansion of $exp\left(x\right)$, $f\left(x\right)$ can be approximated by the PCE with sufficiently large M. In fact, we can calculate the projection $<f\left(x\right),{\varphi }_l\left(x\right)>$ to seek the value of M. When $l>M$, $<f\left(x\right),{\varphi }_l\left(x\right)>\approx 0$. On the other hand, as discussed in Section 4.1, the GPCB can find M automatically by tuning the hyper-parameter $\rho$. Let the experimental design $\mathit{X}$ be the zeros of $H{e}_{10}\left(x\right)$ of Equation (21), i.e., quadrature points corresponding to degree 10; we compare the performance of the GPCB with $\rho$ equal to $0.1,0.45,0.7$, respectively. The results are displayed in Figure 3. Note that the projection value is the absolute value of the true value in the figure for better illustration. It shows that we are able to approximate $f\left(x\right)$ with polynomials up to degree 40. If $\rho =0.45$ for the Mehler kernel, the GPCB almost converges to exact $f\left(x\right)$, whereas $\rho =0.1$ leads to a fast convergence rate and $\rho =0.7$ results in a slow convergence rate.

Figure 3 shows that the $\rho$ has a crucial impact on the performance of the GPCB method, so a tractable method to optimize $\rho$ is needed. A natural criterion is the KL divergence, which can be minimized by finding the optimal hyper-parameters $\rho$. We discussed the KL divergence of the GPCB and the PCE surrogates under the assumption that the PCE surrogate model can approximate the true system to any degree of accuracy. However, the distribution of the real system is usually unknown, which makes the calculation of KL divergence intractable. In fact, it can be easily deduced that the minimization of KL divergence is equivalent to minimizing the negative log marginal likelihood $\Delta$, which (actually is $2\Delta$) reads:

$$\begin{array}{c}\hfill \Delta =Y^T{K}_Y^{-1}Y+log\left[{\left(2\pi \right)}^N\left|K_Y\right|\right]\end{array}$$

(24)

It is important to optimize $\rho$ and ${\sigma }_{ϵ}^2$ to obtain a suitable kernel to get an accurate approximation. Classical methods like gradient-based techniques can be used to search for the optimal $\rho$; however, it may perform poorly because it is locally optimized. As we can see in Equation (23), it is indicated that $\rho$ should take a value between zero and one, so we can propose a global method to solve our optimization problem. The algorithm for generating a GPCB approximation is given in Algorithm 1:

5. Numerical Investigation

In this section, we investigate the GPCB method for various benchmark functions. Firstly, we investigate the same example in Figure 3; however, the noise term is considered, i.e., $y=f\left(x\right)+ϵ=5+x+exp\left(x\right)+ϵ$ is the observation, where $x\sim \mathcal{N}(0,2^2),ϵ\sim \mathcal{N}(0,{0.1}^2)$. Three methods, i.e., GP with the RBF kernel, PCE and GPCB, are implemented. It is necessary to note that the Monte Carlo (MC) sampling strategy is used in the normal GP approaches, and Gaussian quadrature points are introduced in the GPCB approach. The main reason is that the quadrature points are too sparse for the widely-used kernels to capture the local features. For example, in this case, $P=10$ in PCE, then the maximum quadrature point is $10.376$, which is beyond $3{\sigma }_{\mathit{x}}$. In the first set of experiments, let $P=10$ in PCE, which means 11 sample points are used in the experiments; furthermore, 10,000 samples are introduced as test dataset to output the ECDF (empirical cumulative distribution function) and RMSE (root-mean-squared error). The GP algorithm is implemented by the gpml toolbox [38] written in MATLAB with four different kernels, i.e., linear, quadratic, Gaussian and Matérn-3/2 kernels. The comparisons of the results are displayed in Figure 4.

Figure 4a illustrates the point-wise KL divergence between the true value of $f\left(\mathit{x}\right)$ and the predictions on the interval $[-4,4]$ based on Equation (16). It is clear that the distribution of GPCB prediction is statistically closest to the true response, although the GP methods with quadratic, Gaussian and Matérn-3/2 kernels outperform the GPCB at some points. Figure 4b compares the ECDF of y based on the test dataset. It shows that both the PCE and GPCB have a similar ECDF with the true value. Upon closer inspection, which is shown in the magnified subregion, it is obvious that the ECDF of the GPCB is almost exactly the same as the real ECDF, which shows that GPCB has actually captured the feature of $f\left(x\right)$ with high precision. At the same time, the RMSEs of the PCE, linear kernel, quadratic kernel, Gaussian kernel, Matérn-3/2 kernel and GPCB are $1.0570,37.3526,23.6383,25.6044,27.4489,0.4108$, respectively. We have implemented another experiment, which uses the degree P (i.e., the number of experimental design) as the second set of experiments, which are illustrated in Figure 4c. Figure 4c shows that GPCB generally outperforms the ordinary GP with the RBF kernel, which indicates that GPCB performs better with a few (or sparse) training points. It is notable that the PCE and GPCB perform with almost the same precision when the degree is greater than 16. It echoes the idea that PCE and GPCB are statistical equivalent, as we present in Section 4.1.

Similar experiments are conducted with a two-dimensional function, which is expressed as $f\left(\mathit{x}\right)=exp\left(x_1\right)/exp\left(x_2\right)$. Let $x_1,x_2\sim \mathcal{N}(0,1)$, $y=f\left(\mathit{x}\right)+\mathit{ϵ}$ be the real model where $\mathit{ϵ}$ is an independent noise term with a normal distribution $N(0,{0.1}^2{I}_2)$. Unlike the first test function, this test example is a limit state function. Let the maximum degree for each dimension $p_t$ be seven for the PCE method, which makes 64 training points in total. Another dataset of 10,000 independent samples is introduced as the test set to calculate the ECDF and RMSE, as well. Similarly, we have the point-wise KL divergence in the region $[-2,2;-2,2]$ as shown in Figure 5a. It is clear that the GPCB is globally closer to the true distribution than other methods. The GP with quadratic, Gaussian and Matérn-3/2 kernels can approximate the center part well, while the PCE does not seem to perform as well. Figure 5b shows that the six methods except the linear kernel are able to reconstruct the distribution of the prediction, and upon closer observation, we find that the ECDF of the GPCB and Matérn-3/2 kernel are the best approximations among the six methods. We also consider another set of experiments focusing on the number of experimental designs, which equals ${(p_t+1)}^2$ for the 2D function. The RMSEs of the three methods with respect to different $p_t$ are displayed in Figure 5c. It also shows that the PCE and GPCB generally outperform the normal GP approaches, and the GPCB has the best performance.

To summarize, the GPCB generates a surrogate of infinite series, while the PCE can only generate a surrogate with up to $P+1$ polynomials, and they tend to behave with similar precision when P is large enough. A set of sparse quadrature points sampled in the PCE, which are derived from the Gaussian quadrature rule, is a good design for the GPCB. The GPCB with those training points generally performs better than the normal GP methods and PCE. However, the size of such a training set grows dramatically with the dimension ($N={(p_t+1)}^d$ in total), so it is not practical in real-life applications. We aim to present a strategy of sampling from those quadrature points, namely candidate points in the next section, and analyze the performance of our algorithm on the selected points.

5.1. A Random Constructive Design in High Dimensional Problems

As the dimension of a system grows, so do the number of design points of PCE due to a tensor product of quadrature points in each dimension. It is possible that PCE could deal with thousands of points of training data with acceptable computational time; however, it becomes expensive for GP approaches, including our GPCB approach. Monte Carlo sampling techniques can substitute quadrature design; however, these are not always stable. Other sampling strategies like Halton sampling and Latin hypercube sampling [39] are widely used.

In this work, we want to utilize the high accuracy of quadrature points and also want to reduce the massive number of points. Let $\mathit{x}\in {\mathbb{R}}^d$, and $p_t$ is the maximum degree in each dimension, so $p_t+1$ quadrature points are needed in each dimension, which makes the total number of tensor products of quadrature points be $\#\left\{{\mathit{X}}_c\right\}={(p_t+1)}^d$. We seek to find a subset of the candidate design ${\mathit{X}}_c$. Furthermore, we wish to obtain a subset having a good coverage rate in the space. Therefore, we proposed the random definite design in our paper. Note that the LHS design can be extended to a larger interval $(1,{(p_t+1)}^D)$ and can produce points at midpoints (endpoints), so we use the LHS design to sample N indices from the interval. More specifically, we presume that the points in ${\mathit{X}}_c$ are equally important, so we arrange those points with a certain order to get their indices. Then, we sample from the indices with the LHS design, and each index is related to a certain quadrature point. It can be easily implemented by the MATLAB built-in function lhsdesign. The corresponding N points are what we need.

Take a three-dimensional input space as an example, where $x_i\sim \mathcal{N}(0,1),i=1,2,3$. Set $p_t=6$, then $\#\left\{{\mathit{X}}_c\right\}=343$. The candidate design and its subset of 50 points $\mathit{X}$ are illustrated in Figure 6. We can see from the figure that our sampling is sparse in the whole set of candidate points, and it behaves uniformly in dimension one as illustrated in Figure 6b, with similar conclusions in the other two dimensions. When projecting our sampling from dimension three to get Figure 6c, we can see that the selected points almost cover every point of ${\mathit{X}}_c$, which means it has all features in dimensions one and two, i.e., quadrature point values of the two dimensions. It shows that such a method can generate a sparse subset meanwhile guaranteeing the coverage rate in the whole candidate design. We name it the random constructive design.

Now, we want to find out whether these samples retain their capability of accuracy. Firstly, we use the PCE method to test those samples. We will look into the benchmark Ishigami function [40]: $f\left(\mathit{x}\right)=sin\left(x_1\right)+7{sin}^2\left(x_2\right)+0.1x_3^{4}sin\left(x_1\right)$, where four different sampling strategies are compared here. Set $p_t=15$, and the candidate design ${\mathit{X}}_c$ has a size of 4096. The error term $\mathit{ϵ}$ is eliminated in this simulation for the accuracy test. The RMSE are computed on 10,000 independently-sampled data, and the results are presented below in Figure 7. It can be seen that our sample always performs better than other samples. When the number of samples surpasses 900, the RMSE becomes $1.0605\times {10}^{-5}$, which equals the RMSE with the whole candidate points. Therefore, we only select $20\%$ of ${\mathit{X}}_c$ and get the same precision. Furthermore, if we set our precision to be ${10}^{-2}$, only 400 points are needed. This shows that the quadrature points have high precision in numerical calculation. In other words, the points in the candidate set are good points.

Then, the random constructive design is used with the PCE, GP and GPCB methods for the Ishigami function, with the noise term added in the observations. We take the RMSE as a criterion to compare their performance, and the results are illustrated in Figure 8. Figure 8a shows that the GPCB is always better than the PCE method, and they tend to behave the same. However, as the number of sampling points grows, the GP with the quadratic kernel, Gaussian kernel and Matérn-3/2 kernel generally outperform the other methods. We can see that the Ishigami function is a bounded function; therefore, it is likely to fill the whole observation space as the number of samples increases, hence improving the accuracy of the GP method. We plot the ECDF with respect to the three methods when $N=1000$ in Figure 8b. It is clear that the GP with the quadratic kernel, Gaussian kernel and Matérn-3/2 kernel can almost recover the true distribution of the response, which is beyond the capability of the PCE and GPCB.

A six-dimensional problem is being tested with the G-function [41], which is not like the Ishigami function and is unbounded in the domain ${[-\infty ,\infty ]}^6$:

$$\begin{array}{c}\hfill f\left(\mathit{x}\right)=\prod _{i=1}^6\frac{|4x_i-2|+a_i}{1+a_i},where{a}_i=\frac{i-2}{2},\forall i=1,\dots ,6\end{array}$$

(25)

The experiment is performed with the same approaches, and the results are shown below in Figure 9. Figure 9a shows that the GPCB outperforms the PCE and GP with the Gaussian Kernel, and it has similar precision with the quadratic and Matérn-3/2 kernels. We notice that the GPCB is more stable than the PCE method, which behaves badly especially when $N=150,300$. Figure 9b shows that none of these three methods can reconstruct the probability of y very well; however, we can note that the GPCB is still comparatively the closest.

Finally, we are going to present a more complicated model with 15 dimensional functions with the following form:

$$\begin{array}{c}\hfill f\left(\mathit{x}\right)={\mathit{a}}_1^T\mathit{x}+{\mathit{a}}_2^{T}sin\left(\mathit{x}\right)+{\mathit{a}}_3^{T}cos\left(\mathit{x}\right)+{\mathit{x}}^{T}M\mathit{x}.\end{array}$$

(26)

The distribution of $\mathit{x}$ is the product of 15 independent distributions, i.e., $x_i\sim \mathcal{N}(0,1),i=1,\dots ,15$. This function is introduced by the work of O’Hagan, where ${\mathit{a}}_1,{\mathit{a}}_2,{\mathit{a}}_3,M$ are defined in [42]. We can see that this function is dominated by the linear and quadratic term, so it may be well approximated by the low-order PCE model. Let $p_t=3$ in the PCE model; we can see from Figure 10a that the GP with the quadratic kernel performs best among the six methods, while the PCE performs better than the GPCB and other GP methods. On the other hand, the GPCB is always generally better than the GP method except with the quadratic kernel for this function. When $N=1000$, the PCE can generate y, which follows the real distribution according to the ECDF in Figure 10b.

6. Conclusions

This paper has examined two different surrogates of computational models, i.e., polynomial chaos expansion and Gaussian process regression. First, we present a brief review of these two approaches. Next, we discuss the relationship between PCE and GP and find that PCE and GP surrogates are embedded in two isomorphic RKHS. Mercer’s theorem is introduced to generate a kernel based on a PCE basis, by which a new approach is proposed, which we name GPCB. An example shows that with the same experimental design, GPCB tends to retain useful information in a suitable subspace of the RKHS by changing the hyper-parameters, whereas PCE simply sets the information of the residual to zero. We further investigate the approximation performance on two test functions in 1D and 3D, respectively, and their approximation properties are illustrated. In order to deal with the high dimensional scenario, a random constructive design from the quadrature points is used to generate an experimental design. The results give us several directions for choosing models: basically, the GPCB outperforms the PCE, but when the original model can be well approximated by low-order PCE (Figure 10), it seems cumbersome to introduce the GPCB and GP; when the response function is bounded (Figure 8), if we have enough training resources, the GP can be a better choice; when the objective function is unbounded (Figure 4 and Figure 9) or cannot be approximated by finite polynomials (Figure 5), we should probably choose the GPCB.

Future work can extend the family of the Mercer kernel or equivalent kernel (other than the Mehler kernel presented in this paper) beyond a classical approximation method. We can also analyze the experimental design for GP regression in many ways. Although we can see that our sampling method behaves fair enough in the experiments, there is also the opportunity to discover further suitable experimental design schemes to fit different computational purposes, which would be of great interest. The stability of our method will be investigated in future work, i.e., how many points are needed to train a good surrogate and whether our method always produces a suitable design. Furthermore, we can establish closer connections between numerical analysis and statistics via such combinations.