## 2 Kriging Model  
Kriging is widely used to build surrogate models in applied mathematics, statistics, and machine learning. The introduction of its concise theory is given below. For \(d\)-dimensional design data, \(\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{m}]^{T}\,\mathbf{X}\in\Re^{m\times d}\) and objective response \(\mathbf{Y}=[y_{1},\ldots,y_{m}]^{T},\mathbf{Y}\in\Re^{m\times 1}\), the Kriging model is expressed as follows:  
\[y(\mathbf{x})=\mu(\mathbf{x})+z(\mathbf{x}) \tag{1}\]  
The formula consists of two parts. The first part \(\mu(\mathbf{x})\) is the trend function (also called the expectation function). When \(\mu(\mathbf{x})\) is equal to 0, \(\mu\), and \(\sum_{i=1}^{p}\beta f_{i}(\mathbf{x})\), the models are called simple Kriging, ordinary Kriging, and standard Kriging models, respectively, where the function \(f_{i}(x)\) is the regression model, and the parameter \(\beta_{i}\) is the coefficient of the \(i\)-th regression model. Standard Kriging uses the regression function to determine the change in the process mean [38]. The second part in Formula (1) is a stochastic process model established by observing the data and quantifying the correlation of the data. It is often considered to be the realization of the stochastic process \(Z(x)\) with a mean of 0 and a variance of \(\sigma^{2}\). The covariance between design points is shown in Formula (2):  
\[Cov\big{(}Z(x_{i})Z\big{(}x_{j}\big{)}\big{)}=k\big{(}x_{i},x_{j}\big{)}=\sigma ^{2}R\big{(}\mathbf{\theta},\mathbf{x}_{i},\mathbf{x}_{j}\big{)} \tag{2}\]  
where process variance \(\sigma^{2}\) is a scalar and \(\mathbf{\theta}\) is a relative parameter (or called hyper-parameter) of the kernel function. By optimizing \(\mathbf{\theta}\), the correlation between the design points can be adjusted adaptively. The expression of the spatial correlation function \(R\big{(}\mathbf{\theta},\mathbf{x}_{i},\mathbf{x}_{j}\big{)}\) is shown in Equation (3):  
\[R\big{(}\mathbf{\theta},\mathbf{x}_{i},\mathbf{x}_{j}\big{)}=\prod_{k=1}^{d}R_{k}\Big{(} \mathbf{\theta}_{k},\mathbf{x}_{i}^{k}-\mathbf{x}_{j}^{k}\Big{)} \tag{3}\]  
There are many options for spatial correlation functions. However, the most widely used correlation function is the Gaussian function. The expression of the Gaussian model is shown in Equation (4):  
\[R_{k}\Big{(}\mathbf{\theta}_{k},\mathbf{x}_{i}^{k}-\mathbf{x}_{j}^{k}\Big{)}=\exp\left(- \mathbf{\theta}_{k}\Big{|}\mathbf{x}_{i}^{k}-\mathbf{x}_{j}^{k}\Big{|}^{2}\right) \tag{4}\]  
Under the above assumptions, the best linear unbiased estimate of \(y(\mathbf{x})\) is:  
\[\hat{y}(\mathbf{x})=\hat{f}^{T}(\mathbf{x})\hat{\mathbf{\beta}}+\mathbf{r}^{T}(\mathbf{x})\mathbf{R}^{ -1}\Big{(}\mathbf{Y}-\hat{\mathbf{F}}\hat{\mathbf{\beta}}\Big{)} \tag{5}\]  
Parameter \(\mathbf{r}\) and \(\mathbf{R}\) in Formula (5) are:  
\[\mathbf{r}(\mathbf{x})=\{R(\mathbf{x},\mathbf{x}_{1}),R(\mathbf{x},\mathbf{x}_{2}),\cdots,R(\mathbf{x},\bm {x}_{m})\} \tag{6}\]  
Among them, the least square estimation \(\hat{\mathbf{\beta}}\) of \(\mathbf{\beta}\) can be calculated using Equation (7):  
\[\hat{\mathbf{\beta}}=\Big{(}\mathbf{F}^{T}\mathbf{R}^{-1}\mathbf{F}\Big{)}^{-1}\mathbf{F}^{T}\mathbf{ R}^{-1}\mathbf{Y} \tag{7}\]  
In addition, the predicted mean square error \(\hat{s}(\mathbf{x})\) at any point \(\mathbf{x}\) can be obtained by calculating the formula in Equation (8):  
\[\hat{s}(\mathbf{x})=MSE\big{[}\hat{y}(\mathbf{x})\big{]}=\hat{\sigma}^{2}\bigg{\{}1- \Big{[}\hat{f}(\mathbf{x})^{T}\quad\mathbf{r}(\mathbf{x})^{T}\Big{]}\bigg{[} \begin{array}[] {cc}0&\mathbf{F}^{T}\\ \mathbf{F}&\mathbf{R}\end{array} \bigg{]}\bigg{[} \begin{array}{c}\mathbf{f}(\mathbf{x})\\ \mathbf{r}(\mathbf{x})\end{array}