# Step-by-Step Penalized Blind Kriging Methods for Surrogate Modeling

Authors: Ziheng Feng; Min Li; Xiao Wang; Jianxin Zhao; Xinmin Li
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

When building a surrogate model, it is important to identify a proper mean function for the kriging model. The commonly used variable selection method is the penalized blind kriging (PBK) method. But this method could lead to a low time efficiency, which is not suitable for experiments with time-sensitive data. In this paper, we propose three step-by-step approaches for constructing an appropriate mean function to improve the prediction accuracy and time efficiency of the PBK method. Several functions and two engineering examples are used to prove the effectiveness of the proposed methods. From simulation results, we can see that Method 1 (M1) and Method 2 (M2) have been significantly improved in both the prediction accuracy and the time efficiency compared with PBK. Especially, in the Test function, compared with the traditional PBK method, the prediction accuracy of M2 is improved by69.08% and 26.13%, respectively, under the penalty of Lasso and Elastic Net, and the time efficiency of M1 is improved by 85.15% and 90.33%, respectively, under the penalty of Lasso and Elastic Net. In addition, Method 3 (M3) has been significantly improved in prediction accuracy compared with PBK.


---

## Body

## 1. Introduction

The surrogate model has been widely used to replace computation-intensive engineering simulation models or black-box systems [1–3]. There are several commonly used surrogate models, such as kriging, support vector regression, and radial basis functions [4]. Their application fields involve the structural dynamics of aeroengine casings, dynamic reliability estimation of turbine blisk, and computationally expensive constrained optimization problems.The kriging model can be used to fit computer experiments because of its interpolate character. Its error estimation can be extended to many aspects such as failure probability estimation [5], robust design optimization [6], uncertainty analysis [7], and global sequential sampling [8]. As the combination of mean function and stochastic Gaussian process (GP), the kriging model is called ordinary kriging (OK) if using a constant mean to fit the overall trend or universal kriging (UK) if supposing some specified variables in the mean function. However, their prediction accuracy will be reduced if the mean function is specified inappropriately [9]. Computer experiments usually include a large number of input variables, inactive variables which hardly have an impact on the response tending to reduce prediction accuracy; thus, it is a matter to identify active variables into a mean function [10].In this paper, we consider the problem of variable selection for the kriging model. We want to use the method of variable selection to construct an interpretable kriging model which can select the active variables and remove inactive variables with high prediction accuracy and high time efficiency. The commonly used variable selection method in the kriging model is the penalty likelihood method. For the penalty likelihood method, Hung [11] proposed the penalized blind kriging (PBK) approach to select variables. This method, however, suffers from the limitation of low time efficiency and prediction accuracy. In the field of engineering, because of the time sensitive of data, researchers always want to improve time efficiency greatly while maintaining preferable prediction accuracy. Therefore, efficient variable selection methods are called for.The idea here is to improve the time efficiency and prediction accuracy of PBK using step-by-step approaches. Three new approaches are proposed for variable selection. The first method is to select the variables using a linear model at the first and then estimate the parameters in the kriging model. The second method is first to estimate the parameters with the OK model and then select the variables via the penalized kriging model and finally refit the kriging model. The third method is an improvement of PBK, which is a refitting based on the variable selection of PBK. These methods are proved to improve the prediction accuracy and time efficiency by several analytic functions and two examples.The remainder of the paper is organized as follows: in Section2, three new methods are proposed for variable selection. Section 3 illustrates the effectiveness of proposed methods with four simulated examples and compares the prediction accuracy, identification, and time efficiency. Section 4 implements the proposed methods for two engineering examples. Section 5 provides some discussion and concluding remarks.


## 2. Step-by-Step Variable Selection Methods

In this section, we first review some of the details; then, three methods are proposed for variable selection based on PBK methods.


### 2.1. Kriging Preliminaries

The kriging model [12] assumes that the response yx is a realization from the stochastic process:(1)yx=∑i=1pβifix+zx=fxTβ+zx,where fx=f1x,…,fpxT is the basis functions, and β=β1,…,βpT is a vector of unknown coefficients, respectively. The stochastic process zx has mean 0 and covariance σ2R , where Rθ,xi,xj is a correlation function between xi and xj depending on parameter θ. In most engineering applications, the most commonly used correlation function is the stationary Gaussian correlation function Rijθ=Rθ,xi,xj=∏k=1dexp−θkxik−xjk2, which is chosen in our research.Given some sampled pointsxi∈Rd and the responses yi∈R,i=1,…,n, the best linear unbiased predictor (BLUP) can be obtained as follows [12]:(2)y^x=fxTβ^+rxTR−1y−Fβ^,where r=Rθ,x1,x,…,Rθ,xn,xT is a correlation between the unsampled point x and the sampled points xi,i=1,…,n; the matrix F=fx1,…,fxnT is constructed by fx at the sampled points.The log-likelihood function can be written as follows:(3)Lβ,σ2,θ|y=−n2log2π−12logσ2R−12σ2y−FβTR−1y−Fβ.Then, the parameters β, σ2, and θ could be estimated via maximum likelihood estimation (MLE). For given correlation parameters θ, the MLE of β is as follows:(4)β^=FTR−1F−1FTR−1y,and the MLE of σ2 is as follows:(5)σ^2=1ny−Fβ^TR−1y−Fβ^.Since there is no analytical solution for the parametersθ, optimization methods are needed to obtain as follows:(6)θ^=argminθnlogσ^2+logR.Finally, the estimated parametersβ^,σ^2,θ^ are substituted into equation (2) to obtain the predictor of responses for the unsampled points.


### 2.2. Step-by-Step Variable Selection Methods

Although the kriging model has a character of interpolator due to its stochastic process, penalty methods still could be used to select the important variables from the candidates to improve prediction accuracy [13].Two penalty functions are used for variable selection in the PBK model by Hung [11], i.e., Lasso [14] and adaptive Lasso [15]. It estimates the regression coefficients β by minimizing the negative penalized log-likelihood function.(7)Qβ,θ,σ2=y−FβTR−1y−Fβ2σ2+∑j=1pPλ|βj,where Pλβj is the regularization term called “penalty function,” and λ is the regularization parameter. In order to efficiently obtain the estimation of parameters in the penalty likelihood function, Hung [11] proposed the iteratively reweighted least angle regression (IRLARS) algorithm for Lasso penalty and adaptive Lasso penalty.Besides above penalty methods, Elastic Net penalty proposed by Zou and Hastie [16] also could be used for variable selection. It has regularization parameters λ and α. The Elastic Net estimation of the parameters β can be solved by the following formulation:(8)β^EN=argminβy−FβTR−1y−Fβ2σ2+λαβ1+λ1−αβ22,where β1=∑j=1pβj and β22=∑j=1pβj2. Thus, it can be seen that the penalty function of Elastic Net is a convex combination of ridge regression and Lasso. Particularly, the Elastic Net penalty becomes ridge regression when α=0, and the Elastic Net penalty becomes Lasso penalty when α=1. According to the IRLARS algorithm for Lasso penalty, we describe the IRLARS algorithm for Elastic Net penalty in Algorithm 1.Algorithm 1: (IRLARS algorithm for Elastic Net penalty).
Step 1: set the initial values forβ0,θ0,σ20Step 2: withRθ^k−1/σ^2k=CTC,y′=Cy, and F′=CF, solve penalized likelihood function for λ and α obtained by cross validation,β^k+1=argminβy′−F′β2+λα∑j=1pβj+λ1−α∑j=1pβj2Step 3: estimateσ2,θ via maximizing the log-likelihood function (3),σ^2k+1=1/ny−Fβ^k+1TR−1y−Fβ^k+1,θ^k+1=argminθnlogσ^2k+1+logRStep 4: if the convergence is attained, get the estimationβ^,θ^,σ^2; otherwise, return to step 2Regularization parameters are obtained based on 10-fold cross validation, and the estimation of parametersβ,θ,σ2 are obtained based on the IRLARS algorithm in this research.Although PBK effectively reduces the prediction error, further improvements are needed for low time efficiency. Thus, we propose three step-by-step variable selection methods to improve the prediction accuracy and time efficiency. We describe them in detail as follows:


#### 2.2.1. Method One (M1)

We integrate penalty into the linear modelyx=fxTβ+ε firstly, where ε∼N0,σ2; then, we select active variables by solving the penalized likelihood function for β:(9)β^=argminβy−Fβ2+∑j=1pPλβj,where λ is obtained based on cross validation. Then, selected basis functions fi′x,i=1,…,m corresponding nonzero regression coefficients β are substituted into the Kriging model (1) as new variables, and correlation parameters θ are updated as θ′. β′,σ′2,θ′ are estimated by formulas (4)–(6), respectively. Finally, the estimated parameters β^′,σ^′2,θ^′ are used to predict the responses of unsampled points by substituting into formula (2), and the predicted model can be expressed as follows:(10)y^′x=f′xTβ^′+r^′xTR^′−1y−F′β^′.The proposed method M1 is given in Algorithm2.Algorithm 2: (variable selection algorithm for M1).
Step 1: Screen variables using formula (9).Step 2: Select the new basis functionfi′x as new variables. Update the correlation parameter to θ′.Step 3: Estimate the parametersβ′,σ′2,θ′ via formulas (4)–(6).Step 4: The estimated parametersβ^′,σ^′2,θ^′ are used to predict the responses by substituting into formula (10).


#### 2.2.2. Method Two (M2)

Firstly, we estimate the parametersθ^OK and σ^OK2 with the ordinary kriging model yx=μ+zx by solving the following formulas:(11)μ^OK=1nTR−11n−11nTR−1y,(12)σ^OK2=1ny−1nμ^OKTR−1y−1nμ^OK,(13)θ^OK=argminθnlogσ^OK2+logR,where 1n is an n×1 vector of ones. Then, we can select active variables by minimizing the negative penalized log-likelihood function Qβ,θ^OK,σ^OK2 and update correlation parameters. Finally, the new kriging model is fitted similarly to that in M1.The proposed method M2 is given in Algorithm 3.Algorithm 3: (variable selection algorithm for M2).
Step 1: giveμ0, θ0, and σ20 as initial variablesStep 2: use formulas (12) and (13) for optimal model parameters θ^OK,σ^OK2Step 3: withRθ^OK−1/σ^OK2=CTC,y∗=Cy, and F∗=CF solve the penalized likelihood function for λ, β^=argminβy∗−F∗β2+∑j=1pPλβj, where λ is obtained based on cross validationStep 4: the remaining steps of fitting are the same as Steps 2–4 in Algorithm2 and will not be repeated here


#### 2.2.3. Method Three (M3)

Firstly, the PBK method is used to select active variablesβ based on the kriging model via minimizing the negative penalized log-likelihood function (7). It considers correlation between sample points, which is different from M1. Then, the selected variables based on PBK were substituted into the model (1) as new variables, and correlation parameters are updated. Finally, the estimated parameters β^′,σ^′2,θ^′ estimated by formulas (4)–(6) are used to predict the responses of unsampled points by substituting into (2).The variable selection part of M3 is using IRLARS of PBK in Algorithm1. The fitting steps are the same as Steps 2–4 in Algorithm 2 after selecting variables. We will not go into details here.Theoretically, M1 and M2 do not need the iterative step of parameter estimation ofβ such as PBK in variable selection; thus, they can improve time efficiency greatly. In addition, the last step of the proposed three methods is to update θ and refit the kriging model, so the prediction accuracy could be improved.


## 2.1. Kriging Preliminaries

The kriging model [12] assumes that the response yx is a realization from the stochastic process:(1)yx=∑i=1pβifix+zx=fxTβ+zx,where fx=f1x,…,fpxT is the basis functions, and β=β1,…,βpT is a vector of unknown coefficients, respectively. The stochastic process zx has mean 0 and covariance σ2R , where Rθ,xi,xj is a correlation function between xi and xj depending on parameter θ. In most engineering applications, the most commonly used correlation function is the stationary Gaussian correlation function Rijθ=Rθ,xi,xj=∏k=1dexp−θkxik−xjk2, which is chosen in our research.Given some sampled pointsxi∈Rd and the responses yi∈R,i=1,…,n, the best linear unbiased predictor (BLUP) can be obtained as follows [12]:(2)y^x=fxTβ^+rxTR−1y−Fβ^,where r=Rθ,x1,x,…,Rθ,xn,xT is a correlation between the unsampled point x and the sampled points xi,i=1,…,n; the matrix F=fx1,…,fxnT is constructed by fx at the sampled points.The log-likelihood function can be written as follows:(3)Lβ,σ2,θ|y=−n2log2π−12logσ2R−12σ2y−FβTR−1y−Fβ.Then, the parameters β, σ2, and θ could be estimated via maximum likelihood estimation (MLE). For given correlation parameters θ, the MLE of β is as follows:(4)β^=FTR−1F−1FTR−1y,and the MLE of σ2 is as follows:(5)σ^2=1ny−Fβ^TR−1y−Fβ^.Since there is no analytical solution for the parametersθ, optimization methods are needed to obtain as follows:(6)θ^=argminθnlogσ^2+logR.Finally, the estimated parametersβ^,σ^2,θ^ are substituted into equation (2) to obtain the predictor of responses for the unsampled points.


## 2.2. Step-by-Step Variable Selection Methods

Although the kriging model has a character of interpolator due to its stochastic process, penalty methods still could be used to select the important variables from the candidates to improve prediction accuracy [13].Two penalty functions are used for variable selection in the PBK model by Hung [11], i.e., Lasso [14] and adaptive Lasso [15]. It estimates the regression coefficients β by minimizing the negative penalized log-likelihood function.(7)Qβ,θ,σ2=y−FβTR−1y−Fβ2σ2+∑j=1pPλ|βj,where Pλβj is the regularization term called “penalty function,” and λ is the regularization parameter. In order to efficiently obtain the estimation of parameters in the penalty likelihood function, Hung [11] proposed the iteratively reweighted least angle regression (IRLARS) algorithm for Lasso penalty and adaptive Lasso penalty.Besides above penalty methods, Elastic Net penalty proposed by Zou and Hastie [16] also could be used for variable selection. It has regularization parameters λ and α. The Elastic Net estimation of the parameters β can be solved by the following formulation:(8)β^EN=argminβy−FβTR−1y−Fβ2σ2+λαβ1+λ1−αβ22,where β1=∑j=1pβj and β22=∑j=1pβj2. Thus, it can be seen that the penalty function of Elastic Net is a convex combination of ridge regression and Lasso. Particularly, the Elastic Net penalty becomes ridge regression when α=0, and the Elastic Net penalty becomes Lasso penalty when α=1. According to the IRLARS algorithm for Lasso penalty, we describe the IRLARS algorithm for Elastic Net penalty in Algorithm 1.Algorithm 1: (IRLARS algorithm for Elastic Net penalty).
Step 1: set the initial values forβ0,θ0,σ20Step 2: withRθ^k−1/σ^2k=CTC,y′=Cy, and F′=CF, solve penalized likelihood function for λ and α obtained by cross validation,β^k+1=argminβy′−F′β2+λα∑j=1pβj+λ1−α∑j=1pβj2Step 3: estimateσ2,θ via maximizing the log-likelihood function (3),σ^2k+1=1/ny−Fβ^k+1TR−1y−Fβ^k+1,θ^k+1=argminθnlogσ^2k+1+logRStep 4: if the convergence is attained, get the estimationβ^,θ^,σ^2; otherwise, return to step 2Regularization parameters are obtained based on 10-fold cross validation, and the estimation of parametersβ,θ,σ2 are obtained based on the IRLARS algorithm in this research.Although PBK effectively reduces the prediction error, further improvements are needed for low time efficiency. Thus, we propose three step-by-step variable selection methods to improve the prediction accuracy and time efficiency. We describe them in detail as follows:


### 2.2.1. Method One (M1)

We integrate penalty into the linear modelyx=fxTβ+ε firstly, where ε∼N0,σ2; then, we select active variables by solving the penalized likelihood function for β:(9)β^=argminβy−Fβ2+∑j=1pPλβj,where λ is obtained based on cross validation. Then, selected basis functions fi′x,i=1,…,m corresponding nonzero regression coefficients β are substituted into the Kriging model (1) as new variables, and correlation parameters θ are updated as θ′. β′,σ′2,θ′ are estimated by formulas (4)–(6), respectively. Finally, the estimated parameters β^′,σ^′2,θ^′ are used to predict the responses of unsampled points by substituting into formula (2), and the predicted model can be expressed as follows:(10)y^′x=f′xTβ^′+r^′xTR^′−1y−F′β^′.The proposed method M1 is given in Algorithm2.Algorithm 2: (variable selection algorithm for M1).
Step 1: Screen variables using formula (9).Step 2: Select the new basis functionfi′x as new variables. Update the correlation parameter to θ′.Step 3: Estimate the parametersβ′,σ′2,θ′ via formulas (4)–(6).Step 4: The estimated parametersβ^′,σ^′2,θ^′ are used to predict the responses by substituting into formula (10).


### 2.2.2. Method Two (M2)

Firstly, we estimate the parametersθ^OK and σ^OK2 with the ordinary kriging model yx=μ+zx by solving the following formulas:(11)μ^OK=1nTR−11n−11nTR−1y,(12)σ^OK2=1ny−1nμ^OKTR−1y−1nμ^OK,(13)θ^OK=argminθnlogσ^OK2+logR,where 1n is an n×1 vector of ones. Then, we can select active variables by minimizing the negative penalized log-likelihood function Qβ,θ^OK,σ^OK2 and update correlation parameters. Finally, the new kriging model is fitted similarly to that in M1.The proposed method M2 is given in Algorithm 3.Algorithm 3: (variable selection algorithm for M2).
Step 1: giveμ0, θ0, and σ20 as initial variablesStep 2: use formulas (12) and (13) for optimal model parameters θ^OK,σ^OK2Step 3: withRθ^OK−1/σ^OK2=CTC,y∗=Cy, and F∗=CF solve the penalized likelihood function for λ, β^=argminβy∗−F∗β2+∑j=1pPλβj, where λ is obtained based on cross validationStep 4: the remaining steps of fitting are the same as Steps 2–4 in Algorithm2 and will not be repeated here


### 2.2.3. Method Three (M3)

Firstly, the PBK method is used to select active variablesβ based on the kriging model via minimizing the negative penalized log-likelihood function (7). It considers correlation between sample points, which is different from M1. Then, the selected variables based on PBK were substituted into the model (1) as new variables, and correlation parameters are updated. Finally, the estimated parameters β^′,σ^′2,θ^′ estimated by formulas (4)–(6) are used to predict the responses of unsampled points by substituting into (2).The variable selection part of M3 is using IRLARS of PBK in Algorithm1. The fitting steps are the same as Steps 2–4 in Algorithm 2 after selecting variables. We will not go into details here.Theoretically, M1 and M2 do not need the iterative step of parameter estimation ofβ such as PBK in variable selection; thus, they can improve time efficiency greatly. In addition, the last step of the proposed three methods is to update θ and refit the kriging model, so the prediction accuracy could be improved.


## 2.2.1. Method One (M1)

We integrate penalty into the linear modelyx=fxTβ+ε firstly, where ε∼N0,σ2; then, we select active variables by solving the penalized likelihood function for β:(9)β^=argminβy−Fβ2+∑j=1pPλβj,where λ is obtained based on cross validation. Then, selected basis functions fi′x,i=1,…,m corresponding nonzero regression coefficients β are substituted into the Kriging model (1) as new variables, and correlation parameters θ are updated as θ′. β′,σ′2,θ′ are estimated by formulas (4)–(6), respectively. Finally, the estimated parameters β^′,σ^′2,θ^′ are used to predict the responses of unsampled points by substituting into formula (2), and the predicted model can be expressed as follows:(10)y^′x=f′xTβ^′+r^′xTR^′−1y−F′β^′.The proposed method M1 is given in Algorithm2.Algorithm 2: (variable selection algorithm for M1).
Step 1: Screen variables using formula (9).Step 2: Select the new basis functionfi′x as new variables. Update the correlation parameter to θ′.Step 3: Estimate the parametersβ′,σ′2,θ′ via formulas (4)–(6).Step 4: The estimated parametersβ^′,σ^′2,θ^′ are used to predict the responses by substituting into formula (10).


## 2.2.2. Method Two (M2)

Firstly, we estimate the parametersθ^OK and σ^OK2 with the ordinary kriging model yx=μ+zx by solving the following formulas:(11)μ^OK=1nTR−11n−11nTR−1y,(12)σ^OK2=1ny−1nμ^OKTR−1y−1nμ^OK,(13)θ^OK=argminθnlogσ^OK2+logR,where 1n is an n×1 vector of ones. Then, we can select active variables by minimizing the negative penalized log-likelihood function Qβ,θ^OK,σ^OK2 and update correlation parameters. Finally, the new kriging model is fitted similarly to that in M1.The proposed method M2 is given in Algorithm 3.Algorithm 3: (variable selection algorithm for M2).
Step 1: giveμ0, θ0, and σ20 as initial variablesStep 2: use formulas (12) and (13) for optimal model parameters θ^OK,σ^OK2Step 3: withRθ^OK−1/σ^OK2=CTC,y∗=Cy, and F∗=CF solve the penalized likelihood function for λ, β^=argminβy∗−F∗β2+∑j=1pPλβj, where λ is obtained based on cross validationStep 4: the remaining steps of fitting are the same as Steps 2–4 in Algorithm2 and will not be repeated here


## 2.2.3. Method Three (M3)

Firstly, the PBK method is used to select active variablesβ based on the kriging model via minimizing the negative penalized log-likelihood function (7). It considers correlation between sample points, which is different from M1. Then, the selected variables based on PBK were substituted into the model (1) as new variables, and correlation parameters are updated. Finally, the estimated parameters β^′,σ^′2,θ^′ estimated by formulas (4)–(6) are used to predict the responses of unsampled points by substituting into (2).The variable selection part of M3 is using IRLARS of PBK in Algorithm1. The fitting steps are the same as Steps 2–4 in Algorithm 2 after selecting variables. We will not go into details here.Theoretically, M1 and M2 do not need the iterative step of parameter estimation ofβ such as PBK in variable selection; thus, they can improve time efficiency greatly. In addition, the last step of the proposed three methods is to update θ and refit the kriging model, so the prediction accuracy could be improved.


## 3. Simulation

The effectiveness of the proposed methods is proved in terms of simulation performance. We validate the simulation performance of the methods compared with back propagation (BP) neural network, generalized regression (GR) neural network, and PBK method according to time efficiency and prediction accuracy. We use the following indicators to measure time and prediction accuracy.(a)
The root mean square prediction error (RMSPE) is defined as(18)RMSPE=1N∑i=1Nyxitest−y^xitest2,whereN is the number of testing samples. The mean of RMSPE (MRMSPE) of the simulations is calculated to evaluate the prediction performance.(b)
The time of consumption (CPU) shows the uptime of the system.Especially, we need identification to measure their accuracy of variable selection in linear functions using the following indicators:(c)
The average number of inactive effect identified rate (IEIR).(d)
The average number of active effect identified rate (AEIR).(e)
The average size of identified mean function (MEAN).Obviously, when AEIR is larger, MRMSPE, IEIR, and CPU are smaller, and if MEAN is closer to the target, the model is better.Four analytical functions are displayed in the following equations:(1)
A known function [11]:(14)yx=0.4x1+0.3x2+0.2x3+0.1x4+0.05x5+0.01x6+zx.(2)
Test function [17]:(15)yx=9+15∑i=16xi+252∑i=15xixi+1+10∑i=16xi2.(3)
Gramacy and Lee function [13]:(16)yx=expsin0.9x1+0.4810+x2x3+x4.(4)
Borehole function [18]:(17)yx=2πx3x4−x6×logx2x11+2x3x4logx2/x1x12x8−1.A known function (KF) includes twelve-dimensional inputs, where the first six variablesx1,x2,…,x6 decrease the effects on the responses, the remaining variables x7,x8,…,x12, are unrelated (i.e., zero coefficients) to the outputs, and zx has mean 0 and covariance σ2Rθ, where σ2=0.05 and θ=16. The Test function (TF) includes twelve-dimensional inputs, where the last six variables are irrelevant to the outputs. The Gramacy and Lee function (GL) includes six-dimensional inputs, where the last two variables are irrelevant to the outputs. The Borehole function (BF) includes eight-dimensional inputs. And the nonlinearity gradually gets more obvious from the function KF to function BF. In function KF, we consider 60, 80, and 100 sample points for training in order to compare the effect of different sample sizes. In function TF, GL, and BF, we consider 100 sample points for training in order to compare the effect of nonlinearity. 1,000 sample points are randomly selected as testing data. We evaluate the accuracy of prediction, identification, and time efficiency based on 500 simulations.In order to compare the effectiveness of the proposed methods, simulation studies are conducted to evaluate the performance of the BP neural network and GR neural network methods as other surrogate models firstly. The number of neurons in the hidden layer of the BP neural network is the optimal number of neurons among the 15 candidate neurons. The transfer functions from the input layer to the hidden layer and the hidden layer to the output layer choose “logsig” and “purelin,” respectively. The simulation results for four functions are given in Table1.Table 1
Data simulation results of neural network.
FunctionSample SizeBPGRMRMSPECPU(s)MRMSPECPU(s)KFn = 600.110 63008.310.114 52142.63n = 800.108 53509.580.110 52236.98n = 1000.107 73794.200.108 92246.46TFn = 1003.972 75112.7914.469 42429.14GLn = 1000.668 54354.560.672 72306.72BFn = 1004.341 13597.0922.362 22405.54Due to the lack of interpretability of a neural network, neural network methods cannot select variables. In this paper, we will focus on variable selection for the kriging model. We carry out simulation studies using PBK, three methods in Section2 with Lasso penalty and Elastic Net penalty. The first-order polynomial basis function and the stationary Gaussian correlation function are used. And Latin Hypercube designs [19] are used to select sample points for training data in the design domain because they can be produced with minimal computational cost, and the design space is fulfilled well. The simulation results for four functions are given in Tables 2–6, respectively.Table 2
Data simulation results of the function KF with Lasso.
Sample SizeMethodMRMSPEIEIR(%)AEIR(%)MEANCPU(s)n = 60PBK0.067 837.6092.477.802815.65M10.067 555.9396.539.15589.02M20.067 44.3779.205.01687.74M30.066 937.6792.437.812879.55n = 80PBK0.066 720.5391.706.733101.15M10.066 558.1797.009.31618.83M20.066 35.6781.205.21953.33M30.065 720.7791.606.743384.24n = 100PBK0.065 812.5393.176.343525.52M10.065 758.9097.509.38734.71M20.065 07.7782.405.411128.58M30.064 914.2093.176.443881.17Table 3
Data simulation results of the function KF with Elastic Net.
Sample SizeMethodMRMSPEIEIR(%)AEIR(%)MEANCPU(s)n = 60PBK0.067 844.2093.408.2610 315.95M10.067 563.0097.009.601126.76M20.067 06.3381.575.271315.17M30.067 243.4093.638.2210 980.71n = 80PBK0.066 725.1792.277.0510 787.02M10.066 563.1097.279.621277.98M20.066 18.3383.275.501685.08M30.065 825.3792.537.0711 161.00n = 100PBK0.065 717.3793.636.6611 135.97M10.065 663.5797.539.671503.19M20.065 110.7383.405.651914.47M30.065 017.9093.906.7111 761.19Table 4
Data simulation results of the function TF.
MethodLassoElastic NetMRMSPECPU(s)MRMSPECPU(s)PBK0.352 83896.310.355 913 416.90M10.205 1578.450.268 51297.63M20.109 1888.380.262 92286.61M30.108 33924.610.147 614 977.26Table 5
Data simulation results of the function GL.
MethodLassoElastic NetMRMSPECPU(s)MRMSPECPU(s)PBK0.755 52015.580.755 87034.73M10.641 6242.030.640 2766.81M20.633 8600.290.632 61237.27M30.634 22181.550.634 68150.91Table 6
Data simulation results of the function BF.
MethodLassoElastic NetMRMSPECPU(s)MRMSPECPU(s)PBK0.877 92198.330.878 710 522.03M10.710 4367.710.871 31234.14M20.687 0696.510.827 81459.26M30.616 42775.030.796 412 315.72The principal results of these simulations are as follows:(1)
According to Tables2–6, the prediction accuracy of the proposed three methods is better than that of the PBK, BP, and GR neural networks, where M3 is the better one. In addition, their high prediction accuracy is obvious with escalating nonlinearity. Particularly, we take Table 4 as an example; the RMSPE for the predictors with Lasso penalty and Elastic Net in M3 gets low to69.30%=0.3528−0.1083/0.3528and 58.53% smaller than the PBK predictors, respectively.(2)
According to simulation results of the function KF (i.e., Tables2 and 3), the method M1 could identify more variables, which lead to its highest AEIR and IEIR. Although AEIR of M2 is not high, its IEIR is the lowest and MEAN is close to the target. As to M3, its identification is basically the same as PBK. It indicates that M2 and M3 have better identification in variable selection.(3)
In the aspect of time, the efficiency of M1 and M2 is obviously better than that of PBK, BP, and GR neural networks, where M1 is the most efficient. But the efficiency of M3 is little worse than that of PBK. Particularly, we taken=100 in Tables 2 and 3 as an example the time of M1 with Lasso penalty and Elastic Net penalty gets low to 79.16% and 86.50% smaller than PBK, respectively. For other cases, the optimization of time efficiency of M1 is almost higher than 80%.In summary, the BP neural network is mostly worse than the traditional PBK method in the accuracy of prediction, and it is not faster than the PBK method with Lasso penalty in the time efficiency. Although the GR neural network is faster than the PBK method in the time efficiency, it is mostly worse than the traditional PBK method in the accuracy of prediction. We taken=100 in function KF as an example; the RMSPE for BP neural network, the GR neural network and PBK is 0.107 7, 0.108 9, and 0.065 8, respectively, and their time is 3794.20s, 2246.46s, and 3525.52s. In addition, due to the lack of interpretability of BP and GR neural networks, these methods do not have identification. Thus, the neural network is not the preferred candidate method in variable selection methods compared with PBK. Although M1 greatly improves time efficiency by reducing the iterative steps of parameter estimation of β, it identified more inactive variables than PBK, which lead to its high AEIR and IEIR, as the correlation between sample points is not considered in the variable selection. As to M2, it not only maintains the advantage of M1 in time efficiency but also makes up for the shortcomings of M1 by considering the correlation between sample points. Therefore, M2 is superior to M1 in both identification and prediction accuracy and can be used to replace PBK. Although M3 is refitted on the basis of the variable selection of PBK, which leads to its low time efficiency, it has a significant improvement in the prediction accuracy. Therefore, we can use M3 as an improved method of PBK.


## 4. Empirical Application

### 4.1. Circuit Simulation

The proposed methods are applied to the circuit-simulation code [12]. Six input variables x1,…,x6 and one response y are in this experiment. The simulation experiments were conducted 32 times. This data are also used to analyze the relationship between output and input [20]. It can be useful for illustrating the prediction accuracy and time efficiency. The candidate variable set includes all linear effects. To compare the prediction accuracy and time efficiency of different methods, the dataset was separated into two sections: half of the data are used for training and the other half for testing. The experimental results with Lasso and Elastic Net penalty are listed in Tables 7 and 8, respectively.Table 7
Experimental results for circuit simulation with Lasso.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK0.286 9—2.26—M10.220 723.070.3783.63M20.220 723.070.4380.97M30.152 047.022.29-1.33Table 8
Experimental results for circuit simulation with Elastic Net.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK0.274 0—8.65—M10.220 719.450.8190.64M20.220 719.451.0887.51M30.159 141.939.01-4.16From the second and third columns of Tables7 and 8, we could find the prediction accuracy of M1 (0.220 7 and 0.2207), M2 (0.2207 and 0.2207), and M3 (0.1520 and 0.1591) which can achieve better prediction than PBK (0.286 9 and 0.274 0), and the prediction of M1 and M2 are all improved by 23.07% and 19.45%, those of M3 are improved by 47.02% and 41.93% relative to PBK.As displayed in the last two columns of Tables7 and 8, we could find the operation time of M1 (0.37 s and 0.81 s) and M2 (0.43 s and 1.08 s) which are less than those of PBK (2.26 s and 8.65 s) and M3 (2.29 s and 9.01 s), and the approximation efficiency of M1 is improved by 83.63% and 90.64%, those of M2 are improved by 80.97% and 87.51% relative to PBK.Therefore, M1 and M2 can be used as variable selection methods to replacing the traditional PBK method according to prediction accuracy and time efficiency for fitting small samples. And M3 can be used as an improved method of PBK due to better prediction accuracy.


### 4.2. Piston Clap Noise

Piston secondary motion would cause an engine noise which is unwanted. It is the departure of a piston from the nominal motion prescribed by the mechanism. Six variables were used to minimize the piston clap noise. The variables are the cylinder linerx1, location of peak pressure x2, skirt length x3, skirt profile x4, skirt ovality x5, and pin offset x6, respectively.The relevant dataset of the example is derived from Fang et al. [2]. The data include 100 observations, 6 input variables, and the candidate variable set includes all linear effects. In this paper, we divide the 100 data into 5 sets, select one of the sets as the test set for calculating RMSPE, and others as the train sets; this process is repeated 5 times and train sets are different every time. PBK, M1, M2, and M3are compared, and the running time is recorded. The simulation results with Lasso and Elastic Net penalty are shown in Tables 9 and 10, respectively.Table 9
Experimental results for piston clap noise with Lasso.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK1.0069—18.50—M10.831617.412.1188.59M20.546645.715.2071.89M30.635736.8721.66-17.08Table 10
Experimental results for piston clap noise with Elastic Net.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK0.9942—189.45—M10.710828.515.7596.96M20.546645.0216.8991.08M30.546645.02207.76-9.66From the second and third columns of Tables9 and 10, we could find the prediction accuracy of M2 (0.5466 and 0.5466) which can achieve better prediction than PBK (1.0069 and 0.9942), M1 (0.8316 and 0.7108), and M3 (0.6357 and 0.5466), and the prediction of M1 is improved by 17.41% and 28.51%, those of M2 are improved by 45.71% and 45.02%, and those of M3 are improved by 36.87% and 45.02% relative to PBK.As displayed in the last two columns of Tables9 and 10, we can find that the operation time of the M1 (2.11 s and 5.75 s) and M2 (5.20 s and 16.89 s) which are less than those of PBK (18.50 s and 189.45 s) and M3 (21.66 s and 207.76 s), and the approximation efficiency of M1 is improved by 88.59% and 96.96%, and those of M2 are improved by 71.89% and 91.08% relative to PBK.Therefore, M1 and M2 can achieve the goal which replace the PBK method both in terms of prediction accuracy and time efficiency. And M3 can be viewed as an improvement of PBK.


## 4.1. Circuit Simulation

The proposed methods are applied to the circuit-simulation code [12]. Six input variables x1,…,x6 and one response y are in this experiment. The simulation experiments were conducted 32 times. This data are also used to analyze the relationship between output and input [20]. It can be useful for illustrating the prediction accuracy and time efficiency. The candidate variable set includes all linear effects. To compare the prediction accuracy and time efficiency of different methods, the dataset was separated into two sections: half of the data are used for training and the other half for testing. The experimental results with Lasso and Elastic Net penalty are listed in Tables 7 and 8, respectively.Table 7
Experimental results for circuit simulation with Lasso.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK0.286 9—2.26—M10.220 723.070.3783.63M20.220 723.070.4380.97M30.152 047.022.29-1.33Table 8
Experimental results for circuit simulation with Elastic Net.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK0.274 0—8.65—M10.220 719.450.8190.64M20.220 719.451.0887.51M30.159 141.939.01-4.16From the second and third columns of Tables7 and 8, we could find the prediction accuracy of M1 (0.220 7 and 0.2207), M2 (0.2207 and 0.2207), and M3 (0.1520 and 0.1591) which can achieve better prediction than PBK (0.286 9 and 0.274 0), and the prediction of M1 and M2 are all improved by 23.07% and 19.45%, those of M3 are improved by 47.02% and 41.93% relative to PBK.As displayed in the last two columns of Tables7 and 8, we could find the operation time of M1 (0.37 s and 0.81 s) and M2 (0.43 s and 1.08 s) which are less than those of PBK (2.26 s and 8.65 s) and M3 (2.29 s and 9.01 s), and the approximation efficiency of M1 is improved by 83.63% and 90.64%, those of M2 are improved by 80.97% and 87.51% relative to PBK.Therefore, M1 and M2 can be used as variable selection methods to replacing the traditional PBK method according to prediction accuracy and time efficiency for fitting small samples. And M3 can be used as an improved method of PBK due to better prediction accuracy.


## 4.2. Piston Clap Noise

Piston secondary motion would cause an engine noise which is unwanted. It is the departure of a piston from the nominal motion prescribed by the mechanism. Six variables were used to minimize the piston clap noise. The variables are the cylinder linerx1, location of peak pressure x2, skirt length x3, skirt profile x4, skirt ovality x5, and pin offset x6, respectively.The relevant dataset of the example is derived from Fang et al. [2]. The data include 100 observations, 6 input variables, and the candidate variable set includes all linear effects. In this paper, we divide the 100 data into 5 sets, select one of the sets as the test set for calculating RMSPE, and others as the train sets; this process is repeated 5 times and train sets are different every time. PBK, M1, M2, and M3are compared, and the running time is recorded. The simulation results with Lasso and Elastic Net penalty are shown in Tables 9 and 10, respectively.Table 9
Experimental results for piston clap noise with Lasso.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK1.0069—18.50—M10.831617.412.1188.59M20.546645.715.2071.89M30.635736.8721.66-17.08Table 10
Experimental results for piston clap noise with Elastic Net.
MethodRMSPEImproved (%)CPU(s)Improved (%)PBK0.9942—189.45—M10.710828.515.7596.96M20.546645.0216.8991.08M30.546645.02207.76-9.66From the second and third columns of Tables9 and 10, we could find the prediction accuracy of M2 (0.5466 and 0.5466) which can achieve better prediction than PBK (1.0069 and 0.9942), M1 (0.8316 and 0.7108), and M3 (0.6357 and 0.5466), and the prediction of M1 is improved by 17.41% and 28.51%, those of M2 are improved by 45.71% and 45.02%, and those of M3 are improved by 36.87% and 45.02% relative to PBK.As displayed in the last two columns of Tables9 and 10, we can find that the operation time of the M1 (2.11 s and 5.75 s) and M2 (5.20 s and 16.89 s) which are less than those of PBK (18.50 s and 189.45 s) and M3 (21.66 s and 207.76 s), and the approximation efficiency of M1 is improved by 88.59% and 96.96%, and those of M2 are improved by 71.89% and 91.08% relative to PBK.Therefore, M1 and M2 can achieve the goal which replace the PBK method both in terms of prediction accuracy and time efficiency. And M3 can be viewed as an improvement of PBK.


## 5. Conclusion and Discussion

In the field of engineering, we always use a simple surrogate model to fit the complex true model for time saving and interpretability [21]. In order to improve the accuracy of prediction and save the cost of time for the PBK method, three new methods are proposed to obtain the trend function to improve the prediction accuracy and time efficiency. M1 is to select the variables using the linear model at the first and then to estimate the parameters in the kriging model. M2 is first to estimate the parameters with the OK model and then to select the variables via the penalized kriging model and finally refit kriging model. M3 is an improvement of PBK, which is a refitting based on the variable selection of PBK.Several analysis functions with escalating nonlinearity and two engineering examples with different sample sizes are used to test our proposed step-by-step variable selection methods. The experimental results show that M2 which is superior to PBK, BP, and GR neural networks and other methods can be used as a favorable method of variable selection in both prediction accuracy and time efficiency. If we want to pursue higher time efficiency while maintaining certain prediction accuracy, M1 might be a good choice. As well, M3 could be used as an improved method of PBK due to good prediction accuracy.



---

*Source: 4486311-2022-01-21.xml*