计算机试验下Kriging模型选择的比较

Comparison of Model Selection for Kriging Model in Computer Experiments

摘要: Kriging模型是计算机试验的一种常用模型，因具有良好的非线性拟合能力而被广泛使用。本文在一般Kriging模型相关研究的基础上，研究了Kriging模型的变量选择问题，并给出了Elastic Net变量选择方法。与Lasso和adaptive Lasso相比，数值模拟表明Elastic Net变量选择方法能够提高拟合模型的准确性和稳定性。

Abstract: Kriging model is a common model of computer experiment, which is widely used because of its good nonlinear fitting ability. Based on the research of Kriging model, this paper studies the variable selection of universal Kriging model, and gives the variable selection method of Elastic Net. Compared with Lasso and adaptive Lasso, numerical simulation shows that Elastic Net variable selection method can improve the accuracy and stability of the fitting model.

1. 引言

随着科技的不断进步与发展，计算机试验作为物理实验的替代和辅助正在变得越来越流行。计算机试验的一个主要目标是构建一个廉价的元模型。Kriging模型最初是由南非的地质学家Krige在地质统计学中提出并发展的，Sacks等 [1] 于1989年将Kriging模型引入到计算机试验中。从此，Kriging模型作为一种重要的元模型，在农业、工业、化学、生物等领域中有着广泛的应用。

Kriging模型包含平稳高斯过程和均值函数两部分。当均值函数只有一个常数时，称Ordinary Kriging模型；大部分的Kriging模型被假设均值函数部分已存在一些已知变量，这样的模型称为Universal Kriging模型(简称UK)。当模型中包含大量的输入变量时，识别那些(相对较少的)对响应有显著影响的变量是很重要的。Welch等人 [2] 提出可以通过计算机试验进行变量选择，并且指出在保持均值函数不变的情况下，选择对高斯过程有显著影响的变量是非常有意义的。

计算机试验中常用的变量选择方法为贝叶斯方法和惩罚似然方法。对于贝叶斯变量选择，Linkletter等 [3] 提出了高斯过程的贝叶斯选择方法；Huang等 [4] 提出一种盲Kriging模型下对均值函数的变量进行随机搜索的贝叶斯选择方法。对于惩罚似然方法，Li等 [5] 提出高斯过程的惩罚似然变量选择方法，Hung [6] 提出了对均值函数进行Lasso和adaptive Lasso惩罚的似然方法，Zhang等 [7] 提出了一种均值函数的新的惩罚方法，并从贝叶斯观点出发证明了该方法的有效性。

本文研究均值函数的惩罚似然变量选择方法，在Lasso和adaptive Lasso惩罚似然变量的基础上，提出Elastic Net (弹性网络)惩罚。本文的结构如下：第二部分介绍Kriging模型及其参数估计的算法；第三部分介绍了几种变量选择方法及其算法；第四部分通过数据模拟对Kriging模型的几种变量选择方法进行比较；第五部分实例分析，得出结论。

2. Kriging模型简述

已知n个观察点，记为 $X_n={\left(x_1,x_2,\cdots ,x_n\right)}^{\text{T}}$， $Y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$， $y_i=y\left(x_i\right)$ 是点 $x_i$ 对应的输出，其中 $y_i\in R$， $i=1,2,\cdots ,n$。d为设计变量x的维数。一般Kriging模型可写成如下形式：

$y\left(x\right)={\displaystyle \underset{j=1}{\overset{d}{\sum }}{f}_j\left(x\right){\beta }_j}+z\left(x\right)=f^{\text{T}}\left(x\right)\beta +z\left(x\right)=m\left(x\right)+z\left(x\right)$. (2.1)

其中 $m\left(x\right)$ 是均值函数， $f_j\left(x\right),j=0,1,\cdots ,d$ 是已知的基函数， $f\left(x\right)={\left(f_1\left(x\right),\cdots ,f_d\left(x\right)\right)}^{\text{T}}$ 是基向量，d表示回归函数中基函数的个数， ${\beta }_j$ 是未知系数， $\beta ={\left({\beta }_1,\cdots ,{\beta }_d\right)}^{\text{T}}$ 是未知的回归系数向量， $z\left(x\right)$ 是一个高斯随机过程且满足均值为0，即 $E\left(z\left(x\right)\right)=0$，协方差函数为

$\mathrm{cov}\left(z\left(x_i\right),z\left(x_j\right)\right)={\sigma }^{2}c\left(\psi ;x_i,x_j\right),i,j=1,2,\cdots ,n$. (2.2)

在许多计算机试验的文献中 [8]，常用的相关函数为指数相关函数，Matern相关函数及高斯相关函数。本文仅考虑高斯相关函数

$c\left(\psi ;x_i,x_j\right)=\exp \left(-{\displaystyle \underset{i=0}{\overset{d}{\sum }}\frac{\left|x_i-x_j\right|}{{\psi }_i^2}}\right)$. (2.3)

这里的 $\psi$ 是自相关系数， $\psi =\left({\psi }_1,\cdots ,{\psi }_d\right)$。

Santner等 [8] 表明最大似然估计优于交叉验证估计，所以本文使用最大似然估计来估计参数，对数似然函数为

$l\left(\beta ,{\sigma }^2,\psi \right)=-\frac{n}{2}\lg \left(2\pi \right)-\frac{n}{2}\lg \left({\sigma }^2\right)-\frac{1}{2}\lg \left|C\left(\psi \right)\right|-\frac{1}{2{\sigma }^2}{\left(y-F\beta \right)}^{\text{T}}{C}^{-1}\left(\psi \right)\left(y-F\beta \right)$ (2.4)

这里的 $F={\left(f\left(x_1\right),\cdots ,f\left(x_n\right)\right)}^{\text{T}}$ 是 $n\times d$ 的设计矩阵， $C\left(\psi \right)$ 是高斯随机过程 $z\left(x\right)$ 的相关矩阵，第ij个元素为 $c\left(x_i,x_j\right)$。在(2.4)中， $\left(\beta ,{\sigma }^2,\psi \right)$ 均未知。首先假设 $\psi$ 已知，则 $\left(\beta ,{\sigma }^2\right)$ 的估计可以通过最大化(2.4)得到：

$\hat{\beta }={\left(F^{\text{T}}{C}^{-1}\left(\psi \right)F\right)}^{-1}{F}^{-1}{C}^{-1}\left(\psi \right)y$,(2.5)

${\hat{\sigma }}^2=\frac{1}{n}{\left(y-F\beta \right)}^{\text{T}}{C}^{-1}\left(\psi \right)\left(y-F\beta \right)$.(2.6)

此时，(2.4)的最大值为：

$l_{\max }=-\frac{n}{2}\left(1+\lg \left(2\pi \right)\right)-\frac{1}{2}\left(n\log {\hat{\sigma }}^2+\log \left|C\left(\psi \right)\right|\right)$. (2.7)

之后设定 $\psi$ 未知，通过最大化(2.7)得到 $\psi$ 的估计：

$\hat{\psi }=\arg \underset{\psi }{\min }\left\{n\log {\hat{\sigma }}^2+\log \left|C\left(\psi \right)\right|\right\}$.(2.8)

与其它模型相比，Kriging模型可以同时提供预测位置的预测值和预测误差，假定要预测 $x^{\ast }$ 处的函数值， $\hat{y}\left(x^{\ast }\right)$ 表示该点的一个预测值，那么 $\hat{y}\left(x^{\ast }\right)$ 的最佳线性无偏估计(BLUP)为

$\hat{y}\left(x^{\ast }\right)=f^{\text{T}}\left(x^{\ast }\right)\hat{\beta }+c\left(x^{\ast },D\right)C^{-1}\left(y-F\hat{\beta }\right)$.

其中， $c\left(x^{\ast },D\right)$ 是一个 $1\times n$ 的向量，第i个元素是 $c\left(x^{\ast },x_i\right)$。

3. Kriging模型的变量选择方法和算法

考虑Kriging模型的变量选择时，设候选基变量为 $f\left(x\right)={\left(1,f_1\left(x\right),\cdots ,f_p\left(x\right)\right)}^{\text{T}}$， $\beta ={\left({\beta }_0,{\beta }_1,\cdots ,{\beta }_p\right)}^{\text{T}}$， $Z\left(x\right)$ 是均值为0，相关矩阵为 $C\left(\cdot ,\cdot \right)$ 的高斯过程。为了选出 $F\left(x\right)$ 中的重要变量，参数 $\beta$ 需要通过最大化惩罚下面的似然函数估计得到

$Q\left(\beta ,{\sigma }^2,\psi \right)=l\left(\beta ,{\sigma }^2,\psi \right)-{\displaystyle \underset{i=1}{\overset{t}{\sum }}{P}_{\lambda }\left(\left|{\beta }_i\right|\right)}$.

这里的 $l\left(\beta ,{\sigma }^2,\psi \right)$ 是(2.4)， $P_{\lambda }$ 是惩罚函数。

本节首先给出Hung [6] 提出的Lasso方法和adaptive Lasso方法 [9]，然后提出Elastic Net方法。

3.1. Lasso惩罚似然方法

为了估计惩罚Kriging模型中的参数，Hung [6] 提出了重加权最小角度回归算法(IRLARS)，此算法比较容易实施，因此本文所使用的三种惩罚似然方法都要在此算法的框架下进行。

下面详述Lasso惩罚的IRLARS算法，而adaptive Lasso惩罚在Lasso惩罚算法的基础上稍加修改。Lasso惩罚的IRLARS算法如下：

算法1. Kriging模型Lasso惩罚的IRLARS算法

对于adaptive Lasso惩罚的变量选择方法 [9]，是在Lasso的基础上，引入一个已知的权重向量 $w_j$，则adaptive Lasso惩罚下参数 $\beta$ 的参数估计为：

${\hat{\beta }}_{al}^{\left(n\right)}=\underset{\beta }{\arg \min }{\Vert y-{\displaystyle \underset{j=1}{\overset{p}{\sum }}{x}_j{\beta }_j}\Vert }^2+{\lambda }_n{\displaystyle \underset{j=1}{\overset{p}{\sum }}{\hat{w}}_j\left|{\beta }_j\right|}$. (3.1)

这里的 $\hat{w}=1/{\left|\hat{\beta }\right|}^{\gamma }$， $\gamma >0$。在本文研究中， $\hat{w}$ 中的 $\hat{\beta }$ 使用最小二乘估计，即 $\hat{\beta }={\hat{\beta }}_{ols}$， $\gamma =1$。对应的算法只需对算法1中的Step 2稍加修改：

Step 2*：分解 $C{\left({\hat{\psi }}_{\left(l\right)}\right)}^{-1}/{\hat{\sigma }}_{\left(l\right)}^2=R^{\prime }R$，得到 $y^{*}=Ry$， $X^{*}=RX$， $F^{**}=F^{*}/\nu =RF/\nu$，利用交叉验证得到Lasso惩罚下的最优惩罚参数 $\lambda$，求解Lasso问题，

${\hat{\beta }}_{\left(l+1\right)}=\arg \min {\Vert y^{*}-F^{*}\beta \Vert }^2+\lambda {\displaystyle \underset{j=1}{\overset{p}{\sum }}\left|{\beta }_j\right|}$. (3.2)

迭代 ${\hat{\beta }}_j^{\left(l+1\right)}={\hat{\beta }}_j^{*}/{\nu }_j$，这里的 $\nu =\left|{\hat{\beta }}_{ols}\right|$。

3.2. Elastic Net惩罚似然方法

对于Elastic Net惩罚的变量选择方法 [10]，对于固定的非负 ${\lambda }_1$， ${\lambda }_2$，参数 $\beta$ 的Elastic Net估计可通过下式求解，

${\hat{\beta }}_{Enet}=\underset{\beta }{\arg \min }\left\{{\Vert \tilde{Y}-{\tilde{X}}^{\text{T}}\beta \Vert }^2+{\lambda }_2{\Vert \beta \Vert }^2+{\lambda }_1{\Vert \beta \Vert }_1\right\}$. (3.3)

其中， ${\Vert \beta \Vert }^2={\displaystyle \underset{j=1}{\overset{p}{\sum }}{\beta }_i^2}$， ${\Vert \beta \Vert }_1={\displaystyle \underset{j=1}{\overset{p}{\sum }}\left|{\beta }_i\right|}$。由此可见Elastic Net的罚函数是岭回归与Lasso的凸组合。本文研究中，最优惩罚参数 $\lambda$ 依旧通过交叉验证得到， $\alpha =0.7$。

由Elastic Net方法的定义可以看出，当(3.3)式中 ${\lambda }_2=0$ 时，Elastic Net方法就变成了Lasso方法，因此Elastic Net方法结合岭回归与Lasso的优点，既能达到变量选择目的，又具有很好的群组效应。于是只需对算法1中的Step 2稍加修改：

Step 2**：分解 $C{\left({\hat{\psi }}_{\left(l\right)}\right)}^{-1}/{\hat{\sigma }}_{\left(l\right)}^2=R^{\prime }R$，得到 $y^{*}=Ry$， $X^{*}=RX$， $F^{**}=F^{*}/\nu =RF/\nu$，利用交叉验证得到Elastic Net惩罚下的最优惩罚参数 $\lambda$，求解Elastic Net问题，

${\hat{\beta }}_{\left(l+1\right)}=\arg \min {\Vert y^{*}-F^{*}\beta \Vert }^2+{\lambda }_2{\Vert \beta \Vert }^2+{\lambda }_1{\Vert \beta \Vert }_1$.

4. 数据模拟

本节对函数模型做UK模型拟合(UK)，对UK做Lasso惩罚(LUK)，对UK做adaptive Lasso惩罚(ALUK)，和对UK做Elastic Net惩罚(ENUK)，并对这些方法进行比较，最终从变量识别的准确性和预测精度这两方面进行评价。变量识别的准确性用以下三个指标：积极变量识别率的平均(AEIR)；消极变量识别率的平均(IEIR)；积极变量个数的平均(MEAN)。预测精度用以下两个指标来评估：均方根预测误差(RMSPE)的平均值MRMSPE及标准差(sd(RMSPE))。假设有 $n_0$ 个预测点 $x_1^{*},\cdots ,x_{n_0}^{*}$， $\hat{y}\left(x_j^{*}\right)$ 是 $x_j^{*}$ 的BLUP，则RMSPE的计算公式如下RMSPE：

$\text{RMSPE}=\sqrt{\frac{1}{n_0}{\displaystyle \underset{j=1}{\overset{n_0}{\sum }}{\left(y\left(x_j^{*}\right)-\hat{y}\left(x_j^{*}\right)\right)}^2}}$

显然，AEIR越大越好，IEIR，RMSPE，sd(RMSPE)越小越好，MEAN越接近真模型越好。

4.1. 数据模拟1：已知函数

为了评估模型分别加入惩罚前后在预测精度及变量识别方面的性能，我们考虑一个已知线性函数模型 [4]。这个已知函数模型被定义在12维(p = 12)的输入空间 ${\left[0,1\right]}^{12}$ 上，其模型的前六个变量 $x_1,\cdots ,x_6$ 对计算机试验输出的影响逐渐变小，其余变量 $x_7,\cdots ,x_{12}$ 与输出无关(即零系数)。真模型为：

$y\left(x\right)=0.4x_1+0.3x_2+0.2x_3+0.1x_4+0.05x_5+0.01x_6+\epsilon$. (4.1)

这里的 $\epsilon ~N\left(0,{\sigma }_{\epsilon }^2\right)$， ${\sigma }_{\epsilon }=0.05$。响应值y使用(4.1)独立产生，在Matlab上，通过拉丁超立方抽样生成维数p=12，样本量分别为N=50,80,100的样本D，候选变量C由所有的一阶主效应组成。用计算机随机生成的1000个样本点作为测试集G。基于500次的重复试验，模拟结果在表1中给出。

Table 1. Data simulation results of function model (4.1)

表1. 函数模型(4.1)的数据模拟结果

通过表1数据可以发现，从MRMSPE上分析，UK模型都高于LUK，ALUK和ENUK模型，这说明对线性模型进行变量选择可以一定程度上提高模型的预测精度，但LUK，ALUK和ENUK模型的预测精度相差不大。从识别率方面分析，ENUK的积极变量识别率AEIR和识别个数MEAN均高于LUK和ALUK。但LUK的消极变量识别率IEIR低于LUK和ENUK。对于均值函数为多项式模型，模拟结果说明变量选择能够提高预测精度，并且ENUK方法的积极变量识别的准确性方面有明显优势，但消极变量识别率没有优势。

4.2. 数据模拟2：钻孔函数

考虑钻孔函数 [11] 如下：

$y\left(x\right)=2\pi x_3\left(x_4-x_6\right){\left\{\log \left(\frac{x_2}{x_1}\right)\left(1+2\frac{x_3{x}_4}{\log \left(x_2/x_1\right)x_1^2{x}_8}+\frac{x_3}{x_5}\right)\right\}}^{-1}$ (4.2)

输入空间是一个矩形区间[0.05, 0.015] × [100, 5000] × [63,070, 115,600] × [990, 1110] × [63.1, 116] × [700, 820] × [1120, 1680] × [9855, 12,045]，拟合模型形式为 $\hat{y}\left(x\right)=m\left(x\right)+z\left(x\right)$，设定均值函数形式为 $m\left(x\right)={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_8{x}_8$。在Matlab上，通过拉丁超立方抽样生成n = 100，d = 8的样本D，计算机随机生成1000个样本点作为测试集G，计算均方根预测误差RMSPE。同样地，对(4.2)做UK, LUK, ALUK, ENUK几种方法的拟合。重复进行500次试验，并计算每种情况下的MRMSPE以及标准差sd (RMSPE)得到表2。

Table 2. Data simulation results of the functional model (4.2)

表2. 函数模型(4.2)的数据模拟结果

从表2数据可以看出，从MRMSPE和sd (RMSPE)的数据分析，UK模型的MRMSPE和sd (RMSPE)远高于LUK, ALUK和ENUK模型，这说明对非线性模型进行变量选择也可以一定程度上提高模型的预测精度，但LUK, ALUK和ENUK模型的预测精度相差不大。从变量选择方面分析，ENUK的识别个数MEAN高于LUK和ALUK。对于均值函数为多项式模型，模拟结果表明变量选择能够提高预测精度，并且ENUK方法的变量识别有优势。这说明在真实模型为非线性模型的条件下，变量选择大大提高了参数估计的准确性与稳定性。

5. 实例分析：活塞拍击噪声

本节应用本文提出的对Kriging模型做变量选择的各种方法，对一个活塞拍击噪声的实例进行分析。活塞拍击是由活塞二次运动引起的一种不必要的发动噪声。通过计算机试验，改变六种因素来减少排挤噪声。这六种因素分别为活塞与缸套之间的间隙下x₁，峰值压力的位置x₂，活塞裙部长度x₃，裙部型线x₄，裙部椭圆轮廓x₅，活塞销偏置x₆。实例的相关数据集来源于Huang [4]，数据包括100个观测值，6个输入变量，候选变量集C包括，所有的线性主效应，二次主效应，正交多项式编码下的两因素交叉效应 [12]，因此C一共包含72个基变量。本研究进行5折交叉验证，每次将100个观测值随机选出80个数据作为训练集用于构建模型，剩下的20个作为预测集用于计算RMSPE，比较UK模型，OK模型，以及分别作Lasso，adaptive Lasso，Elastic Net惩罚前后的MRMSPE，sd (RMSPE)，MEAN。比较UK模型，以及分别做Lasso，adaptive Lasso，Elastic Net惩罚后得到的LUK，ALUK，ENUK的MRMSPE，sd(RMSPE)， MEAN。模拟结果如表3所示。

Table 3. Data simulation results of a piston slap noise example

表3. 活塞拍击噪声实例的数据模拟结果

从表3的数据可以看出，ENUK的MRMSPE和sd (RMSPE)远比UK, LUK和ALUK模型小；模型变量的个数远低于UK，而且ENUK模型的变量个数与LUK和ALUK相差不大。通过这个实例结果显示，使用ENUK方法可以有效简化模型的同时还可以降低均方预测误差。

6. 结论

本文在一般Kriging模型相关研究的基础上，研究了Kriging模型的变量选择问题，并给出了Elastic Net变量选择方法。模拟结果和实例验证表明，Elastic Net变量选择方法与Lasso和adaptive Lasso相比，Elastic Net能够提高拟合模型的准确性和稳定性。