青岛大学数学与统计学院，山东 青岛

收稿日期：2023年2月21日；录用日期：2023年3月20日；发布日期：2023年3月27日

摘要

本文研究多任务Kriging模型的变量选择问题，并给出多种稀疏化惩罚下多任务Kriging的变量选择算法。数值模拟及实例分析表明，相比单任务的Kriging变量选择，多任务模式能显著提高计算效率而不失模型拟合的准确性；相比LMC及卷积模型，多任务稀疏化Kriging能有效提取任务间的共性信息，极大节约计算成本同时提高预测精度。

关键词

多响应Kriging模型，模型选择，多任务学习，惩罚函数

Research and Application on Variable Selection in Multi-Task Kriging Model

Jie Ji, Chenchen Zou^*

School of Mathematics and Statistics, Qingdao University, Qingdao Shandong

Received: Feb. 21^st, 2023; accepted: Mar. 20^th, 2023; published: Mar. 27^th, 2023

ABSTRACT

We study the variable selection in multi-task Kriging model and develop the algorithms for commonly used penalizations. In numerical simulations, our multi-task penalized approach achieves higher computational efficiency without loss of accuracy and stability compared to the single-task approach. In real data application, multi-task penalized Kriging effectively captures shared features among tasks and thus reduces computational burden compared with the LMC and CONV models.

Keywords:Multi-Response Kriging Model, Variable Selection, Multi-Task Learning, Penalization

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

计算机试验通过建立相对低成本的元模型来模拟复杂物理问题中的输入输出关系 [1]。Kriging模型作为一种经典的元模型，其变量选择被广泛的研究和应用 [2]。研究者先后提出Lasso Kriging [3]，惩罚盲Kriging [4]，惩罚盲似然Kriging [5] [6]，随机搜索盲Kriging [7] 及稀疏回归Kriging模型 [8] 等。

目前变量选择的研究普遍基于单任务Kriging模型。对于具有多响应的系统，逐个应用单任务式的变量筛选，一方面容易忽略任务间的关联及共享信息，另一方面未能利用多任务学习计算效率上的优势。多任务的Kriging模型受到越来越多的关注，但截至目前，多任务Kriging变量选择的研究相对有限 [9]。

变量选择对于多任务Kriging模型的意义不仅限于核心特征筛选。多任务学习在Kriging趋势函数部分的正则化惩罚可有效提取任务间共享信息，相比常用于处理多任务高斯过程的线性协同区域模型 (LMC) [10] 和卷积模型(CONV) [11]，稀疏化的多任务Kriging能以更低的计算成本拟合任务间的关联性，同时具备相当甚至更好的预测精度。

本文研究多任务Kriging模型趋势函数部分的变量选择方法。内容安排如下：第二部分介绍多任务 Kriging模型；第三部分介绍多任务Kriging模型变量选择算法；第四部分为数值模拟研究；第五部分为实例应用及分析；第六部分为研究总结。

2. 多任务Kriging模型简述

设系统中有m个输出如下：

$Y_t\left(x\right)=f^{\text{T}}\left(x\right){\beta }_t+z_t\left(x\right),t=1,\cdots ,m$ (1)

其中 $f\left(x\right)={\left(f_1\left(x\right),\cdots ,f_p\left(x\right)\right)}^{\text{T}}$ 是已知的基函数， ${\beta }_t={\left({\beta }_{t1},\cdots ,{\beta }_{tp}\right)}^{\text{T}}$ 是待估的系数向量， $z_t\left(x\right)$ 是零均值、相互独立的高斯过程 $\left(t=1,\cdots ,m\right)$ 。对于任意给定的两点 $x$ 和 $x^{\prime }$ ， $Y_t\left(x\right)$ 和 $Y_t\left(x^{\prime }\right)$ 的协方差为

$K_t\left(x,x^{\prime }\right)=cov\left(Y_t\left(x\right),Y_t\left(x^{\prime }\right)\right)=cov\left(z_t\left(x\right),z_t\left(x^{\prime }\right)\right)={\sigma }_t^2{R}_t\left(x,x^{\prime }\right)$ (2)

其中 $R_t\left(x,x^{\prime }\right)=\exp \left(-{\displaystyle {\sum }_{h=1}^d{\theta }_{th}}\left|x_h-x_{h^{\prime }}\right|\right)$ ， ${\theta }_t$ 是自相关系数。本文采取约束极大似然法估计参数，对数似然函数为

$\ln L=-\frac{1}{2}{\displaystyle \underset{t=1}{\overset{m}{\sum }}\left[{\left(Y_t-F\left(X_t\right){\beta }_t\right)}^T{K}_t^{-1}\left(Y_t-F\left(X_t\right){\beta }_t\right)+\ln \left|K_t\right|\right]}$ (3)

其中， $X_t={\left(x_{t1},\cdots ,x_{t{n}_t}\right)}^{\text{T}}$ 表示 $n_t$ 组观测点的输入， $F\left(X_t\right)=\left(f_j{x}_{ti}\right),i=1,\cdots ,n_t;j=1,\cdots ,p$ 表示 $n_t\times p$ 的基函数矩阵， $K_t$ 表示第t个任务 $n_t\times n_t$ 的协方差矩阵， $t=1,\cdots ,m$ 。给定 ${\theta }_t$ ，最大化(3)式得到 $\left({\beta }_t,{\sigma }_t^2\right)$ 的估计

${\hat{\beta }}_t\left({\theta }_t\right)={\left(F^T\left(X_t\right)K_t^{-1}F\left(X_t\right)\right)}^{-1}F{\left(X_t\right)}^{\text{T}}{K}_t^{-1}{Y}_t$ (4)

${\hat{\sigma }}_t^2\left({\theta }_t\right)=\frac{{\left(Y_t-F\left(X_t\right){\hat{\beta }}_t\right)}^{\text{T}}{K}_t^{-1}\left(Y_t-F\left(X_t\right){\hat{\beta }}_t\right)}{n_t}$ (5)

将(4) (5)代入(3)式关于 ${\theta }_t$ 最大化可得

${\hat{\theta }}_t=\arg {\min }_{{\theta }_t}\left\{\left(n_t-p\right)\ln {\hat{\sigma }}_t^2\left({\theta }_t\right)+\ln \left|K_t\right|+\ln \left|F^T\left(X_t\right)K_t^{-1}F\left(X_t\right)\right|\right\}$ (6)

在预测点 $x^{\ast }$ ，最佳线性无偏估计为

${\hat{y}}_t\left(x^{\ast }\right)=f^T\left(x^{\ast }\right){\hat{\beta }}_t+k_{t^{\ast }}^T{K}_t^{-1}\left(Y_t-F\left(X_t\right){\hat{\beta }}_t\right),t=1,\cdots ,m$ (7)

其中， $k_{t^{\text{*}}}=K_t\left(X_t,x^{\text{*}}\right)\in {\text{R}}^{n_t\times 1}$ ， $t=1,\cdots ,m$ 。

3. 多任务Kriging模型变量选择算法

本文通过惩罚似然法估计模型参数。记 $\beta =\left({\beta }_1,\cdots ,{\beta }_m\right)$ ， $P\left(\lambda ;\beta \right)$ 表示由超参数λ调整的惩罚函数，如Multi-task Lasso [12] 、L₂₁-norm [13] 及Dirty Model [14]。多任务Kriging的变量选择算法步骤如下：

算法3.1

Step 1：设置初值 ${\hat{\sigma }}_t^{\left(0\right)}$ ， ${\hat{\theta }}_{th}^{\left(0\right)}$ ， $h=1,\cdots ,d$ ，代入(2)式得到协方差矩阵 $K_t$ ， $t=1,\cdots ,m$ ；

Step 2：Cholesky分解 $K_t^{-1}=C_t^{\text{T}}{C}_t$ ，令 ${\tilde{Y}}_t=C_t{Y}_t$ ， ${\tilde{F}}_t=C_{t}F\left(X_t\right)$ ，最小化下式求解 ${\hat{\beta }}_t^{\left(1\right)}$ ， $t=1,\cdots ,m$ ，

${\displaystyle \underset{t=1}{\overset{m}{\sum }}{\Vert {\tilde{Y}}_t-{\tilde{F}}_t{\beta }_t\Vert }_2^2}+P\left(\lambda ;\beta \right)$ (8)

$\lambda$ 利用交叉验证调节；

Step 3：将 ${\hat{\beta }}_t^{\left(1\right)}$ 代入(3)、(6)式求得更新后的 ${\hat{\sigma }}_t^{\left(1\right)}$ ， ${\hat{\theta }}_t^{\left(1\right)}$ ， $t=1,\cdots ,m$ ；

Step 4：重复Step 2，Step 3直至收敛。

Multi-task Lasso [12] 通过最小化下式得到 ${\hat{\beta }}_t$ ， $t=1,\cdots ,m$ 。

$\underset{{\beta }_t}{\min }{\displaystyle \underset{t=1}{\overset{m}{\sum }}{\Vert {\tilde{Y}}_t-{\tilde{F}}_t{\beta }_t\Vert }_2^2}+\lambda {\displaystyle \underset{t=1}{\overset{m}{\sum }}{\Vert {\beta }_t\Vert }_1}$ (9)

其中， ${\Vert {\beta }_t\Vert }_1={\displaystyle {\sum }_{j=1}^p\left|{\beta }_{tj}\right|}$ ，对应修改Step 2里 $P\left(\lambda ;\beta \right)=\lambda {\displaystyle {\sum }_{t=1}^m{\Vert {\beta }_t\Vert }_1}$ 即可。

L₂₁-norm方法 [13] 对 $\beta$ 逐行稀疏化筛选，为m个任务选择共享特征，同时实现下述目标函数：

$\underset{{\beta }_t}{\min }{\displaystyle \underset{t=1}{\overset{m}{\sum }}{\Vert {\tilde{Y}}_t-{\tilde{F}}_t{\beta }_t\Vert }_2^2}+\lambda {\displaystyle \underset{t=1}{\overset{m}{\sum }}\sqrt{{\displaystyle \underset{j=1}{\overset{p}{\sum }}{\beta }_{tj}^2}}}$ (10)

对应修改Step 2里 $P\left(\lambda ;\beta \right)=\lambda {\displaystyle {\sum }_{t=1}^m\sqrt{{\displaystyle {\sum }_{j=1}^p{\beta }_{tj}^2}}}$ 即可。

Dirty Model方法 [14] ：Dirty Model的关键思想是将 $\beta$ 分解为P和Q两个分量， ${\lambda }_1$ 控制P上的群稀疏正则化， ${\lambda }_2$ 控制Q上的稀疏正则化，鼓励所有任务选择相同的一组特征(通过组稀疏组件)，目标函数表示为：

$\underset{{\beta }_t}{\min }{\displaystyle \underset{t=1}{\overset{m}{\sum }}{\Vert {\tilde{Y}}_t-{\tilde{F}}_t\left(P_t+Q_t\right)\Vert }_2^2}+{\lambda }_1{\Vert P\Vert }_{1,\infty }+{\lambda }_2{\Vert Q\Vert }_1$ (11)

其中 ${\Vert P\Vert }_{1,\infty }=\max {\Vert P_j\Vert }_{1,\infty },{\Vert Q\Vert }_1={\displaystyle {\sum }_{t=1}^m{\displaystyle {\sum }_{j=1}^p\left|Q_{tj}\right|}}$ ，对应修改Step 2里 $P\left(\lambda ;\beta \right)={\lambda }_1{\Vert P\Vert }_{1,\infty }+{\lambda }_2{\Vert Q\Vert }_1$ 。

4. 数值模拟研究

我们对单任务Lasso惩罚(Lasso)，多任务Lasso惩罚(MTL)，L₂₁-norm惩罚(L₂₁)以及Dirty Model惩罚(Dirty)下的Kriging模型，拟合效果进行模拟比较,。软件采用MATLAB，以工具箱MALSAR [15] 为主进行多任务计算。基本原理见第三部分。评估从以下5个指标进行：积极变量识别率均值(AEIR)、消极变量识别率均值(IEIR)、模型长度均值(MEAN)、均方根预测误差平均值(MRMSPE)以及均方根预测误差标准差(sd (RMSPE))。AEIR越大越好；IEIR，MRMSPE，sd (RMSPE)越小越好；MEAN越接近真实值越好。

模拟:本研究根据以下三种模型生成响应数据，分别是：

模型I：

$y_{\text{t}}={\beta }_{t,0}+{\beta }_{t,1}{x}_1+{\beta }_{t,2}{x}_2+{\beta }_{t,3}{x}_3+z_t\left(x\right),t=1,2$ (13)

模型II：

$y_t={\beta }_{t,0}+{\beta }_{t,1}{x}_1+{\beta }_{t,2}{x}_2+{\beta }_{t,3}{x}_3+{\beta }_{t,4}{x}_1^2+{\beta }_{t,5}{x}_2^2+{\beta }_{t,6}{x}_3^2+z_t\left(x\right),t=1,2$ (14)

模型III：

$\begin{array}{l}{y}_t={\beta }_{t,0}+{\beta }_{t,1}{x}_1+{\beta }_{t,2}{x}_2+{\beta }_{t,3}{x}_3+{\beta }_{t,4}{x}_1^2+{\beta }_{t,5}{x}_2^2+{\beta }_{t,6}{x}_3^2+{\beta }_{t,7}{x}_1{x}_2\\ +{\beta }_{t,8}{x}_1{x}_3+{\beta }_{t,9}{x}_2{x}_3+z_t\left(x\right),t=1,2\end{array}$ (15)

其中， $x={\left(x_1,\cdots ,x_8\right)}^{\text{T}}$ ，对于每个模型，回归系数独立地从 $\left[-20,-10\right]\cup \left[10,20\right]$ 随机选择。此外，我们随机选择 $p_k\times {\rho }_k$ ，其中 ${\rho }_p=1/3$ 是一个模拟参数，相应的回归系数从 $\left[-10,0\right]\cup \left[0,10\right]$ 中随机选择。 $z\left(x\right)$ 是一个中心平稳多元高斯过程，其相关函数由exp相关函数

$R_t\left(x,x^{\prime }\right)=\exp \left(-{\displaystyle \underset{h=1}{\overset{d}{\sum }}{\theta }_{th}\left|x_h-x_{h^{\prime }}\right|}\right)$ (16)

给出，参数 $\theta$ 在[0, 2]上随机生成，过程方差被设置为1。

样本从[0, 1]上通过拉丁超立方抽样取得，对于上述设置运行50次，训练样本量设为50，预测样本量设为500，结果见表1。

表1. 模拟多任务Kriging模型的变量选择

单任务Kriging相比多任务Kriging，预测精度上总体相差不大，多任务对复杂情况表现更稳定，多任务比单任务更省时。

多任务之间比较而言，预测精度上MTL对不同复杂情况的预测更加稳定，而Dirty 和L₂₁模型会随着模型复杂程度增加，其预测精度均出现越来越差和较大的波动；变量识别率上，三者对积极变量的识别率都在80%以上，对消极变脸的识别率，MTL明显更有优势，这是由于Dirty和L₂₁都有强制为不同的任务选择共同变量的特性，这就会出现在任务间存在差异时共享信息中包含不适用本任务信息的情况，从而导致它们的IEIR相对较高；计算上，Dirty Model的运算时间远高于MTL、L₂₁。

由上表可知，当真实模型为多响应高斯过程时，多任务Kriging模型变量选择能够在保证预测精度的同时减少运算时间。相较于Dirty和L21，MTL更稳定，且更能好地简化模型。

5. 实例分析

我们将MTL模型、LMC模型以及CONV模型用于在伊朗德黑兰370栋住宅公寓数据集 [16]，进行分析预测，旨在提供一个关于在设计阶段或施工初期估算任何给定城市的新住房价格的代理模型，该模型将显示与新建筑单元销售价格相关的影响因素。每组观测值由27个输入和2个输出组成，其中，输入变量为房地产单元的物理和财务属性(x₁~x₈)以及相关的经济变量(x₉~x₂₇)，输出变量为实际销售价格以及实际建筑费用。

我们采用5折交叉验证调节惩罚函数，超过75%的任务中明显不重要的输入包括市颁发的建筑许可证总建筑面积x₁₃，时间分辨率下的贷款利率x₁₉和私营部门在建筑开始时的平均建筑成本x₂₀，房地产单元的物理和财务属性(x₁~x₈)为显著重要变量，这与Rafiei [16] 的研究结果一致。

我们比较MTL、LMC、CONV三种模型在测试集上的均方根预测误差平均值。LMC由Multi-output-Gaussian-Process [17] 包实现，CONV由multigp [18] 实现，均采用默认设置。MTL的预测效果最好，LMC次之。模型训练时间上，MTL对多个任务的共性提炼相当于模拟高斯过程不同响应间的相关性，相比LMC和CONV对任务间相关性的挖掘方式，大大提高了运算效率，同时预测更加精准(见表2)。

表2. 三个模型的预测结果及训练模型时间比较

6. 总结

在单任务Kriging模型相关研究的基础上，本文研究了多任务Kriging模型的变量选择问题，并给出了该模型的与建筑销售和实际费用实例相关的试验。模拟结果和实例验证表明，与LMC和COMV相比，多任务Kriging模型变量选择方法能够提高拟合模型的准确性并有效降低其运算成本。多任务Kriging模型能够筛选出对新房市场更具影响的因素，并筛除不重要变量，简化模型，充分借鉴多任务之间的相关信息，本研究使用相对较少样本来训练模型，可有效减轻收集数据所带来的成本问题，可操作性更强，同时可训练出泛化性能更好的新住房价格的代理模型。

文章引用

纪 洁,邹晨晨. 多任务Kriging变量选择的研究与应用
Research and Application on Variable Selection in Multi-Task Kriging Model[J]. 应用数学进展, 2023, 12(03): 1224-1230. https://doi.org/10.12677/AAM.2023.123124

参考文献