## 3 Proposed Algorithm
### Nonstationary Kriging
After the input samples are divided by utilizing the criteria function in Equation (13), the correlation matrix based on the given \(N\) samples is constructed as  
\[\mathbf{R}\big{(}\zeta^{i},\zeta^{i}\big{)}=\left[ \begin{array}{ccc}R\big{(} \zeta^{1},\zeta^{1}\big{)}&\ldots&R\big{(}\zeta^{1},\zeta^{N}\big{)}\\ \vdots&\ddots&\vdots\\ R\big{(}\zeta^{N},\zeta^{1}\big{)}&\ldots&R\big{(}\zeta^{N},\zeta^{N}\big{)} \end{array} \right] \tag{14}\]  
The kriging model parameters are computed using the best linear unbiased estimator (BLUE) [22] given by  
\[ \begin{array}{l}\beta=\big{(}\mathbf{F}^{\mathrm{T}}\mathbf{R}^{-1}\mathbf{F }\big{)}^{-1}\mathbf{F}^{\mathrm{T}}\mathbf{R}^{-1}\mathbf{y}\\ \sigma^{2}=\Big{[}\big{(}\mathbf{y}-\mathbf{F}\beta\big{)}^{\mathrm{T}} \mathbf{R}^{-1}(\mathbf{y}-\mathbf{F}\beta)\Big{]}/N\end{array}  \tag{15}\]  
The hyperparameters are obtained by solving an optimization problem that minimizes the concentrated log-likelihood function as follows:  
\[\psi(\theta)=\frac{1}{2}\ln\big{|}\mathbf{R}(\theta)\big{|}+\frac{N}{2}\ln \big{(}\sigma^{2}(\theta)\big{)} \tag{16}\]  
In solving the problem in Equation (16), the pattern search, which is one of the global optimization algorithms, is adopted. After obtaining hyperparameters, the predicted response and its variance at the unsampled input point \(\zeta\) are computed as  
\[ \begin{array}{l}y(\zeta)\simeq\mathbf{f}^{\mathrm{T}}(\zeta)\beta+\mathbf{r }^{\mathrm{T}}(\zeta)\mathbf{R}^{-1}(\mathbf{y}-\mathbf{F}\beta)\\ \hat{\sigma}^{2}(\zeta)=\sigma^{2}\Bigg{(}1-\left[\begin{array}{cc}\mathbf{ f}^{\mathrm{T}}(\zeta)&\mathbf{r}^{\mathrm{T}}(\zeta)\end{array}\right]\left[ \begin{array}{cc}0&\mathbf{F}^{\mathrm{T}}\\ \mathbf{F}&\mathbf{R}\end{array}\right]^{-1}\Bigg{[}\begin{array}{c}\mathbf{ f}(\zeta)\\ \mathbf{r}(\zeta)\end{array}\right]\Bigg{)}\end{array}  \tag{17}\]  
where \(\mathbf{r}(\zeta)\) is the correlation vector between input \(\zeta\) and given \(N\) samples, and its \(i\)-th component is  
\[r_{i}(\zeta) =R\Big{(}\zeta,\zeta^{i}\Big{)}=\prod_{k=1}^{M}\Bigg{(}\frac{ \overline{\hat{l}_{k}}(\zeta)l_{k}\Big{(}\zeta^{i}\Big{)}}{\overline{\hat{l}_{ k}^{2}}(\zeta)+\bar{l}_{k}^{2}\big{(}\zeta^{i}\Big{)}}\Bigg{)}^{1/2}\exp\left(- \frac{\big{(}\zeta-\zeta^{i}\big{)}^{2}}{\overline{\hat{l}_{k}^{2}}(\zeta)+ \bar{l}_{k}^{2}\big{(}\zeta^{i}\big{)}}\right) \tag{18}\]  
where \(\overline{\hat{l}_{k}}(\zeta)\) is the \(k\)-th length-scale for the unsampled input \(\zeta\). Figure 5: Comparison of nonstationary and stationary correlation function. (**a**) Nonstationary case; (**b**) Stationary case. One difference between stationary kriging and the proposed method is that, in the case of stationary kriging, one can estimate the response for any input using Equation (17). However, for nonstationary kriging, the correlation function in Equation (18) cannot be constructed, since the length scale function for unsampled input is unknown. As a result, the response prediction in Equation (17) is not directly applicable. To tackle this issue, support vector machine (SVM) -based classification for unsampled input is applied and will be discussed in the next subsection.