多精度计算机试验的高维代理模型

High-Dimensional Surrogate Modeling for Multi-Fidelity Computer Experiments

摘要: 多精度计算机试验结合不同精度的数据，以降低计算成本并提高预测精度。然而，在高维问题中，传统方法面临计算复杂度高、数据稀疏性强等挑战。本文提出一种基于变换可加高斯(TAG)过程的高维多精度代理模型，并结合嵌套设计，在精度层级间构建层级化误差修正关系。该方法有效提升了高维情形下的建模精度和计算效率。数值实验表明，本文提出的方法在高维问题中优于传统的多精度高斯过程模型，适用于高维情形下的多精度计算机试验的代理建模。

Abstract: Multi-fidelity computer experiments integrate data of varying fidelity levels to reduce computational costs and improve prediction accuracy. However, traditional methods face challenges in high-dimensional problems due to high computational complexity and severe data sparsity. This paper proposes a high-dimensional multi-fidelity surrogate model based on the transformation additive Gaussian (TAG) process and incorporates the nested design to establish a hierarchical error correction relationship. The proposed approach effectively enhances modeling accuracy and computational efficiency in high-dimensional cases. Numerical experiments demonstrate that this method outperforms traditional multi-fidelity Gaussian process models in high-dimensional settings, making it well-suited for surrogate modeling of high-dimensional multi-fidelity computer experiments.

1. 引言

在计算科学与工程领域，计算机试验被广泛用于模拟复杂物理[1]、化学[2]、生物[3]等系统，以辅助科学研究与工程设计。然而，随着计算问题的复杂度和维度的增加，高精度计算模型往往需要大量计算资源，使得直接进行高精度仿真变得极为昂贵甚至不可行。为了解决这一问题，多精度建模(multi-fidelity modeling)逐渐成为一种有效的替代方案，其核心目标是通过融合不同精度数据源实现高效建模与预测。它结合了高精度模型和低精度模型的优势，以较低的计算成本获得可靠的预测结果。其中，多精度高斯过程(Multi-Fidelity Gaussian Process, MFGP)是最常见的建模方法之一。例如，在碳纤维/环氧树脂层合板的冲击损伤预测中，高精度模型采用晶体塑性有限元模拟单晶尺度的滑移系与位错演化，而低精度模 型基于宏观各向异性损伤本构方程。通过高斯过程回归建立微观–宏观应变场的映射关系，仅在关键载荷步调用高精度数据，最终将误差控制在8%以内，显著提升计算效率[4]。多精度高斯过程基于经典的高斯过程(Gaussian Process, GP) [5]回归方法，利用不同精度层级的数据来提高预测精度，同时降低计算成本。Kennedy和O’Hagan [6]提出的经典协方差加权(Co-Kriging)方法假设低精度模型与高精度模型之间满足线性相关关系，并利用贝叶斯框架进行联合建模。该模型在多精度计算机试验的各个领域做出了重大贡献[7]-[9]，为多精度建模奠定了基础。基于此框架，Perdikaris [10]提出了一种结合非线性自回归的多精度高斯过程，能够在多精度数据中学习复杂的非线性和空间相关性。Qian [11] [12]，Le Gratiet [13] [14]和Ji [15]等人研究了高效后验预测和贝叶斯不确定性量化的建模策略。这些方法通常基于以下假设：低精度模型能够捕捉高精度模型的主要趋势，而残差或偏差可以通过统计方法进行修正。与此同时，多精度模型也在不断地发展与应用，在工程优化领域，NASA研究中心于2017年开发出EGO-MF (Efficient Global Optimization with Multi-Fidelity)算法[16]，其核心在于引入了贝叶斯主动学习机制，通过EI (Expected Improvement)函数动态分配高、低精度样本的采样位置。麻省理工学院(MIT)研究团队于2016年提出的Deep MFGP (Deep Multi-Fidelity Gaussian Processes) [17]，通过堆叠式神经网络自动提取输入–输出空间的非线性映射特征，并耦合高斯过程实现不确定性量化。尽管多精度高斯过程在低维问题上取得了成功，但在高维场景下面临诸多挑战。首先，高斯过程的计算复杂度随样本数量和输入维度的增加呈指数级增长，导致计算成本过高。而多精度模型需同时处理高、低精度数据的协方差矩阵，进一步加剧计算负担。此外，高维空间的样本分布稀疏性和非线性交互效应导致传统协方差核(如径向基核)难以建模此类复杂结构。因此，如何在高维多精度计算机试验中构建高效且准确的代理模型仍然是一个重要的研究问题。Tsilifis等人[18]提出了一种利用低维投影的贝叶斯高斯过程建模方法，以降低高维建模的困难。但线性投影假设可能无法捕捉复杂的非线性降维关系，从而导致信息损失。同时基于贝叶斯推断和哈密顿蒙特卡洛的采样方法仍有着较高的计算量。Wu等人[19]提出了一种多精度分层神经过程(Multi-Fidelity Hierarchical Neural Processes, MF-HNP)模型，用于解决高维非嵌套多精度数据的代理建模问题。但模型需要大量数据来提升预测性能，若数据不足或质量差，模型表现可能受限。

本文提出了一种基于变换可加高斯(Transformation Additive Gaussian, TAG)过程的新型多精度建模方法。TAG过程由Lin和Roshan [20]于2020年提出，是一种适用于高维问题的高斯过程建模方法。与传统的高斯过程模型不同，TAG过程通过对响应变量进行适当的变换，使复杂的高维函数在变换后的尺度上近似满足加性结构，从而将高维问题降解为多个低维问题，有效降低了高维优化问题的复杂度。在本文中，我们结合嵌套设计，在不同精度水平的数据之间构建层级化的误差修正关系，使得最终的高精度预测可以通过低精度模型及其修正函数高效获取。该方法不仅能够降低计算成本，还能够在高维问题中保持较好的预测精度。本文的主要贡献如下：提出了一种适用于高维空间的多精度代理建模框架，该框架结合嵌套设计与TAG过程，有效利用不同精度级别的数据，实现精度层级间的信息传递；通过数值实验验证该方法在高维多精度计算机试验中的有效性，并与现有方法进行对比分析。本文的其余部分安排如下：第二部分介绍变换可加高斯过程；第三部分详细描述所提出的高维多精度代理模型的构建；第四部分进行数值实验分析，验证方法的有效性；第五部分总结全文。

2. TAG过程介绍

在许多计算机试验和高维复杂系统的建模问题中，传统的高斯过程方法由于全局协方差矩阵的高维优化、参数识别困难以及对变量间复杂交互效应的处理不足，常常难以满足大规模、高维数据分析的需求。为此，Lin和Roshan [20]提出了一种基于响应变量转换的模型：TAG过程。其核心思想在于对响应变量进行适当变换，使得原本非线性且含有交互效应的函数在转换后的尺度上尽可能呈现出加性结构，从而实现对各自单变量效应的分解，同时通过引入适当的修正项捕捉残留的交互效应。

考虑一个确定性函数

$y=f\left(x\right),$

其中$y$ 为响应变量，$x=\left(x_1,x_2,\cdots ,x_p\right)$ 为$p$ 维输入变量。传统GP方法直接对$y$ 建模，但在无噪声情形下直接拟合高维函数不仅需要大量数据，而且计算复杂度较高。TAG过程首先引入变换函数$g(\cdot )$ ，如Box-Cox变换，其定义为[21]

$g_{\lambda }\left(y\right)=\{\begin{array}{l}\frac{y^{\lambda }-1}{\lambda },\text{if }\lambda \ne 0\hfill \\ \log y,\text{ if }\lambda =0\hfill \end{array}$

其中$\lambda$ 为需要估计的参数。通过这种变换，可以期望使得转换后的响应满足加性，即假设有

$g_{\lambda }\left(y\right)=\mu +{\displaystyle \sum _{k=1}^p{z}_k}\left(x_k\right)+\epsilon \left(x\right),$

其中$\mu$ 为全局均值，$z_k\left(x_k\right)$ 表示仅依赖于第$k$ 个自变量的单变量效应，而$\epsilon \left(x\right)$ 为较小的误差项，满足$\epsilon \left(x\right)~N\left(0,{\sigma }^2\right)$ 。类似于传统的高斯过程，我们可以使用贝叶斯框架得到预测函数。假设$z_k\left(x_k\right)$ 为零均值的高斯过程，则其具有正态先验

其中${\tau }_k^2$ 为方差，$R_k\left(h\right)=Cor\left\{z_k\left(x_k\right),z_k\left(x_k+h\right)\right\}$ 为平稳相关函数，常用的选择是高斯相关函数

$R_k\left(h\right)=\exp \left\{-\frac{h^2}{s_k^2}\right\},$

其中$s_k$ 为尺度参数。定义整体方差为

${\tau }^2={\displaystyle \sum _{k=1}^p{\tau }_k^2},$

并设权重参数

${\omega }_k=\frac{{\tau }_k^2}{{\tau }^2},满足{\displaystyle \sum _{k=1}^p{\omega }_k}=1,$

则整体加性协方差函数可以构造为

$R\left(h\right)={\displaystyle \sum _{k=1}^p{\omega }_k}{R}_k\left(h_k\right).$

因此，TAG过程在转换后的尺度上可表示为

$g_{\lambda }\left(y\right)=\mu +z\left(x\right)+\epsilon \left(x\right),其中z\left(x\right)~GP\left(0,{\tau }^{2}R(\cdot )\right).$

给定数据${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^n$ ，经过变换得到$g_{\lambda }\left(y\right)={\left(g_{\lambda }\left(y_1\right),\cdots ,g_{\lambda }\left(y_n\right)\right)}^{\prime }$ ，定义$\delta ={\sigma }^2/{\tau }^2$ (作为“nugget”参数，反映模型逼近误差)，则$g_{\lambda }\left(y\right)$ 具有分布

且$z\left(x\right)$ 的后验分布为

$z\left(x\right)|y~N\left(\widehat{z}\left(x\right),{\tau }^2\left\{1+\delta -r{\left(x\right)}^{\prime }{\left(R+\delta I\right)}^{-1}r\left(x\right)\right\}\right),$

其中$\widehat{z}\left(x\right)=r{\left(x\right)}^{\prime }{\left(R+\delta I\right)}^{-1}\left(g\left(y\right)-\mu 1\right)$ ，$r\left(x\right)={\left(R\left(x-x_1\right),\cdots ,R\left(x-x_n\right)\right)}^{\prime }$ ，$R={\left(R\left(x_i,x_j\right)\right)}_{n\times n}$ ，$I$ 为单位矩阵，$1$ 表示元素全为1的列向量。

下面通过贝叶斯框架得到未知参数的估计。数据变换$g_{\lambda }\left(y\right)$ 在非信息先验下，对参数$\theta =\left(\mu ,{\tau }^2,\delta ,\omega ,s,\lambda \right)$ 的后验分布与边际似然有关，其对数似然包含项

$-\frac{n}{2}\log {\tau }^2-\frac{1}{2}\log \left|R+\delta I\right|-\frac{1}{2{\tau }^2}{\left[g_{\lambda }\left(y\right)-\mu 1\right]}^{\prime }{\left(R+\delta I\right)}^{-1}\left[g_{\lambda }\left(y\right)-\mu 1\right],$

以及由于变换引入的雅可比校正项

${\displaystyle \prod _{i=1}^n{y}_i^{\lambda -1}}.$

通过最大化此边际似然，可以获得$\mu$ 与${\tau }^2$ 的估计，即

$\widehat{\mu }=\frac{1^{\prime }{\left(R+\delta I\right)}^{-1}{g}_{\lambda }\left(y\right)}{1^{\prime }{\left(R+\delta I\right)}^{-1}1},$

${\widehat{\tau }}^2=\frac{1}{n}{\left[g_{\lambda }\left(y\right)-\widehat{\mu }1\right]}^{\prime }{\left(R+\delta I\right)}^{-1}\left[g_{\lambda }\left(y\right)-\widehat{\mu }1\right],$

而其他参数(如$\delta ,\lambda ,s,\omega$ )则需要通过数值优化求得。得益于模型的加性结构，初始参数可以利用回归拟合(例如借助R包$mgcv$ 中的backfitting算法[22])以及一维GP拟合快速获得，从而将原在$2p+2$ 维空间中的全局优化问题降至低维问题，显著提高了计算效率。

尽管TAG过程经过适当的响应变换后能够较好地近似加性结构，但在实际问题中，转换后的响应仍可能存在部分高阶交互效应。为此，TAAG (Transformed Approximately Additive Gaussian)过程在TAG过程的基础上对误差项$\epsilon \left(x\right)$ 进行进一步建模，假设其服从具有平稳协方差结构的高斯过程，即

$\epsilon \left(x\right)~GP\left(0,{\sigma }^{2}L(\cdot )\right),$

其中$L(\cdot )$ 一般采用乘积型高斯相关函数

$L\left(h\right)={\displaystyle \prod _{k=1}^p\exp }\left\{-\frac{h_k^2}{{\gamma }_k^2}\right\},$

${\gamma }_k$ 为误差过程的尺度参数。这样，TAAG过程整体可以表示为

其中

${\nu }^2={\tau }^2+{\sigma }^2,\eta =\frac{{\sigma }^2}{{\tau }^2+{\sigma }^2}\in \left[0,1\right].$

参数$\eta$ 控制了非加性成分的贡献：当$\eta$ 趋近于0时，TAAG过程退化为TAG过程；而较大的$\eta$ 则表明转换后仍存在显著的交互效应。

在参数估计方面，为避免高维优化的复杂性，通常先利用TAG过程得到$\lambda ,\omega ,s$ 的初步估计，再固定这些参数，对TAAG过程中混合参数$\eta$ 及误差过程尺度参数$\gamma$ 进行一维或二维优化，从而大幅降低整体优化问题的维度。

模型的预测部分依然利用GP的条件分布性质构造。给定观测数据，预测点$x$ 处转换后响应的后验分布为

其中预测均值为

预测方差为

其中$l\left(x\right)=L\left(x-x_1,\cdots ,x-x_j\right)$ 。由于$g_{\lambda }(\cdot )$ 通常为非线性函数，为便于逆变换和不确定性量化，实际中常采用预测中位数作为最终预测值，并通过数值方法计算可信区间。

总之，TAG过程通过对响应变量进行Box-Cox等变换，将原始复杂的高维函数转化为各变量主效应的加性组合，不仅改善了预测精度，还通过参数${\omega }_k$ 直观地反映了各变量的重要性；而TAAG过程在此基础上引入平滑的误差过程，通过参数$\eta$ 捕捉转换后残留的交互效应，进一步提高了模型的灵活性与不确定性量化能力。两者均采用经验贝叶斯方法进行参数估计，并利用加性结构大幅降低全局优化的计算复杂度，使得在大样本与高维数据下仍具有良好的适用性。因此，TAG与TAAG过程不仅在预测性能上优于传统GP模型，同时在解释性、可视化及变量敏感性分析方面也提供了有力支持。

3. 基于TAG过程的多精度建模

3.1. 嵌套设计与样本构造

以两精度计算为例，我们假设存在两个计算机代码，其中低精度代码(LC)运行成本低但精度有限，而高精度代码(HC)精度更高但运行代价较大。基于HC和LC数据构建多精度代理模型的研究受到广泛关注，其基本思路通常是先利用一对嵌套的空间填充设计(如嵌套拉丁超立方设计[23] [24])对HC与LC分别进行采样，再基于LC数据构建初步预测模型，并利用HC数据对模型进行修正和校准[6] [25]-[27]。在此过程中，嵌套设计的思想对多精度建模具有重要意义，因为它能够确保低精度样本包含高精度样本，即满足$\left(D_{HC}\subset D_{LC}\right)$ ，从而在两个层次之间建立稳固的结构关联，降低数据不匹配带来的误差。

嵌套设计的核心在于在同一设计空间内，通过构造一个空间填充设计作为LC样本，并从中选取一部分点作为HC样本来实现。这样，高精度样本不仅在空间上均匀分布，还能完全包含于低精度样本中，有效保证了两层数据之间的一致性，为后续模型校正提供了便利。此外，为了进一步优化样本构造，通常会采用最大熵设计、空间填充度等准则[28] [29]来挑选HC样本，以确保这些点在设计空间中具有良好的代表性和均匀覆盖性。这种嵌套设计结构不仅在一次性设计中具有优势，更适用于序贯设计策略。Xiong等人[30]通过逐步扩充LC样本的同时，保持HC样本作为其子集，可以在每一阶段都评估代理模型的预测性能，并根据交叉验证误差或条件交叉验证误差决定是否需要新增样本。这样的数据驱动式采样过程，不仅节省了HC计算资源，还能动态调整设计规模，使得多精度代理模型既能充分利用廉价的低精度信息，又能通过高精度样本及时校正模型偏差，从而达到精度和稳定性的提升。因此，嵌套设计不仅保证了多精度数据的一致性，还为TAG过程在高维情形下的有效建模提供了充足的样本支持和均匀的空间覆盖，进而提升模型的解释性与预测能力。

3.2. 基于TAG过程的多精度建模

在多精度计算机试验中，我们希望通过低精度数据与少量高精度数据的结合，构建一个能够逼近高精度结果的高效代理模型。传统方法通常采用GP对低精度数据建模，并进一步建模高低精度数据之间的偏差模型。然而，在高维情况下，GP建模的计算复杂度较高，同时难以有效捕捉变量间的复杂非线性关系。本节以TAAG过程为例，构建多精度模型，以提高建模精度并降低计算成本。

首先，以两精度为例，利用嵌套设计生成低精度样本$D_l=\left\{\left(x_1,y_1^{(l)}\right),\cdots ,\left(x_{n_l},y_{n_l}^{(l)}\right)\right\}$ 和高精度样本$D_h=\left\{\left(x_1,y_1^{(h)}\right),\cdots ,\left(x_{n_h},y_{n_h}^{(h)}\right)\right\}$ ，满足$D_h\subset D_l$ ，其中$y^{(l)}=\left(y_1^{(l)},\cdots ,y_{n_l}^{(l)}\right)$ 与$y^{(h)}=\left(y_1^{(h)},\cdots ,y_{n_h}^{(h)}\right)$ 分别表示低精度响应和高精度响应。利用TAAG过程对低精度数据进行建模，即

$g_{{\lambda }^{(l)}}\left(y^{(l)}\left(x\right)\right)={\mu }^{(l)}+{\displaystyle \sum _{k=1}^p{z}_k^{(l)}}\left(x_k\right)+{\epsilon }^{(l)}\left(x\right),$

其中$g_{{\lambda }^{(l)}}$ 表示Box-Cox变换函数，${\lambda }^{(l)}$ 为变换参数。$z_k^{(l)}\left(x_k\right)~\text{GP}\left(0,{\tau }_k^{(l)2}{R}_k^{(l)}(\cdot )\right)$ 包含尺度参数$s_k^{(l)}$ ；残差项${\epsilon }^{(l)}\left(x\right)~\text{GP}\left(0,{\sigma }^{(l)}{}^2{L}^{(l)}(\cdot )\right)$ ，采用乘积高斯相关函数$L^{(l)}\left(h\right)={\displaystyle \prod _{k=1}^p\exp }\left(-h_k^2/{\gamma }_k^{(l)2}\right)$ ，用来捕捉未被可加结构解释的交互效应。参数集${\theta }^{(l)}=\left({\mu }^{(l)},{\nu }^{(l)2},{\lambda }^{(l)},{\omega }^{(l)},s^{(l)},{\tau }^{(l)2},{\gamma }^{(l)},{\eta }^{(l)}\right)$ 通过最大化边际似然函数估计，正如Lin [20]所介绍的那样，可以使用一些策略来简化计算。例如固定${\lambda }^{(l)},{\omega }^{(l)},s^{(l)}$ 为TAG过程中的估计值，构造一个GP来拟合$g_{{\lambda }^{(l)}}\left(y^{(l)}\left(x\right)\right)$ ，获得${\gamma }^{(l)}$ 的估计值，给参数${\eta }^{(l)}$ 一个Beta先验等，使得最后所需要估计的参数仅有${\mu }^{(l)}$ ，${\nu }^{(l)2}$ 和${\eta }^{(l)}$ ，进一步由贝叶斯公式有

${\widehat{\eta }}^{(l)}=\arg {\min }_{{\eta }^{(l)}}\log \left|\left(1-{\eta }^{(l)}\right){\widehat{R}}^{(l)}+{\eta }^{(l)}{L}^{(l)}\right|+n\log {\widehat{\nu }}^{(l)2}-2\log \left\{{\eta }^{(l)}\left(1-{\eta }^{(l)}\right)\right\},$

其中

${\widehat{\mu }}^{(l)}=\frac{1^{\prime }{\left\{\left(1-{\eta }^{(l)}\right){\widehat{R}}^{(l)}+{\eta }^{(l)}{L}^{(l)}\right\}}^{-1}{g}_{{\widehat{\lambda }}^{(l)}}\left(y^{(l)}\right)}{1^{\prime }{\left\{\left(1-{\eta }^{(l)}\right){\widehat{R}}^{(l)}+{\eta }^{(l)}{L}^{(l)}\right\}}^{-1}1},$

${\widehat{\nu }}^{(l)2}=\frac{1}{n}{\left(g_{{\widehat{\lambda }}^{(l)}}\left(y^{(l)}\right)-{\widehat{\mu }}^{(l)}1\right)}^{\prime }{\left\{\left(1-{\eta }^{(l)}\right){\widehat{R}}^{(l)}+{\eta }^{(l)}{L}^{(l)}\right\}}^{-1}\left(g_{{\widehat{\lambda }}^{(l)}}\left(y^{(l)}\right)-{\widehat{\mu }}^{(l)}1\right),$

其中${\widehat{R}}^{(l)}$ 为$R^{(l)}$ 将${\omega }^{(l)}$ 和$s^{(l)}$ 代入后的近似，$L^{(l)}={\left(L^{(l)}\left(x_i-x_j\right)\right)}_{n_l\times n_l}$ 。

变换后低精度函数具有后验分布

其中$\widehat{g°y^{(l)}}\left(x\right)$ ，$V^{(l)}\left(x\right)$ 分别为$\widehat{g°y}\left(x\right)$ 和$V\left(x\right)$ 对应低精度的形式。我们进一步构建偏差修正模型，假设高精度响应可以表示为

$y^{(h)}\left(x\right)=y^{(l)}\left(x\right)+\delta \left(x\right),$

我们将修正项利用TAAG过程进行建模，由于经典的Box-Cox变换仅适用于响应变量为正的情况，这里采用其拓展形式[21]

$g_{\lambda }\left(y\right)=\{\begin{array}{l}\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1},\text{if }{\lambda }_1\ne 0\hfill \\ \log \left(y+{\lambda }_2\right),\text{ if }{\lambda }_1=0\hfill \end{array}$

从而

$g_{{\lambda }_{\delta }}\left(\delta \left(x\right)\right)={\mu }_{\delta }+{\displaystyle \sum _{k=1}^p{z}_{\delta ,k}}\left(x_k\right)+{ϵ}_{\delta }\left(x\right),$

其中${\lambda }_{\delta }=\left({\lambda }_{\delta 1},{\lambda }_{\delta 2}\right)$ 为未知参数，在实际操作中，我们需要选择${\lambda }_{\delta 2}$ ，满足$\delta +{\lambda }_{\delta 2}>0$ 。其参数利用高低精度残差$\left\{y^{(h)}\left(x_i\right)-y^{(l)}\left(x_i\right)\right\}$ 来优化，具有分布

其中$\widehat{g°y_{\delta }}\left(x\right)$ ，$V_{\delta }\left(x\right)$ 分别为$\widehat{g°y}\left(x\right)$ 和$V\left(x\right)$ 对应偏差$\delta$ 的形式。最终的高精度预测模型为

$\widehat{y}\left(x\right)={\widehat{y}}^{(l)}\left(x\right)+\widehat{\delta }\left(x\right).$

这一建模策略充分结合了低精度数据的计算优势和高精度数据的精度优势，在保证计算成本可控的同时，实现高维数据下的高精度预测。实验结果表明，该方法在多个高维测试问题上均取得了优于传统高斯过程建模的效果。

此外，我们可以基于分布构造高精度预测模型的可信区间。从转换后的分布$g_{{\lambda }^{(l)}}°y^{(l)}\left(x\right)|y^{(l)}$ 和$g_{{\lambda }_{\delta }}°y_{\delta }\left(x\right)$ 中抽取独立样本，通过逆变换函数(如Box-Cox逆变换$g_{\lambda }^{-1}(\cdot )$ )将样本映射至原始响应尺度，得到低精度预测样本$y_i^{(l)}=g_{\lambda }^{-1}\left(g_{\lambda }\left(y_i^{(l)}\right)\right)$ 和修正项样本${\delta }_i^{(l)}=g_{{\lambda }_{\delta }}^{-1}\left(g_{{\lambda }_{\delta }}\left({\delta }_i^{(l)}\right)\right)$ ，将二者逐点叠加生成高精度预测样本$y_i^{(h)}=y_i^{(l)}+{\delta }_i^{(l)}$ ，从而$100\left(1-\alpha \right)\%$ 可信区间为

$\left[y_{[(\alpha /2)n]}^{(h)},y_{[(1-\alpha /2)n]}^{(h)}\right],$

其中为$y_i^{(h)}$ 的次序统计量中第i个值，$[\cdot ]$ 表示取整函数。

我们可以把这一过程拓展到多精度试验中。假设存在$L$ 个精度层级，其中第$l$ 层级($l=1,\cdots ,L$ )的响应函数为$y^{(l)}\left(x\right)$ ，且$y^{(l)}\left(x\right)$ 为计算成本最高的高精度模型。采用嵌套设计生成层级化试验点集${\left\{D^{(l)}\right\}}_{l=1}^L$ ，满足严格的嵌套性：低层级设计$D^{(l)}$ 完全包含于高层级设计$D^{(l+1)}$ 中，即$D^{(1)}\subset D^{(2)}\subset \cdots \subset D^{(L)}$ 。每个层级$l$ 包含$n_l$ 个样本点($n_1<n_2<\cdots <n_L$ )。针对各层级数据，首先建立单层TAAG过程。对于第$l$ 层级的响应$y^{(l)}\left(x\right)$ ，通过Box-Cox变换$g_{{\lambda }^{(l)}}(\cdot )$ 将其映射至可加空间：

$g_{{\lambda }^{(l)}}\left(y^{(l)}\left(x\right)\right)={\mu }^{(l)}+{\displaystyle \sum _{k=1}^p{z}_k^{(l)}}\left(x_k\right)+{ϵ}^{(l)}\left(x\right),$

类似地，$z_k^{(l)}\left(x_k\right)~\text{GP}\left(0,{\tau }_k^{(l)2}{R}_k^{(l)}\left(x_k,{x^{\prime }}_k;s_k^{(l)}\right)\right)$ 为独立的高斯过程分量，${\epsilon }^{(l)}\left(x\right)~\text{GP}\left(0,{\sigma }^{(l)}{}^2{L}^{(l)}(\cdot )\right)$ ，参数集${\theta }^{(l)}=\left({\lambda }^{(l)},{\tau }^{(l)2},{\sigma }^{(l)2},{\omega }^{(l)},s^{(l)},{\gamma }^{(l)},{\eta }^{(l)}\right)$ 仍通过最大化边际似然函数估计。其中，可以将初始值通过嵌套设计传递，高层级的参数${\theta }^{(l+1)}$ 初始化为低层级的估计值${\theta }^{(l)}$ ，以充分利用低层级数据的先验信息。为融合多层级信息，进一步构建层级间误差修正模型。假设相邻层级间存在线性映射关系，定义第$l+1$ 层级响应为：

$y^{(l+1)}\left(x\right)=y^{(l)}\left(x\right)+{\delta }^{(l)}\left(x\right),$

修正项${\delta }^{(l)}\left(x\right)$ 采用TAAG模型建模：

$g_{{\lambda }_{\delta }^{(l)}}\left({\delta }^{(l)}\left(x\right)\right)={\mu }_{\delta }^{(l)}+{\displaystyle \sum _{k=1}^p{z}_{\delta ,k}^{(l)}}\left(x_k\right)+{ϵ}_{\delta }^{(l)}\left(x\right),$

其参数${\theta }_{\delta }^{(l)}$ 基于高层级残差$\left\{y^{(l+1)}\left(x_i\right)-{\widehat{y}}^{(l)}\left(x_i\right)\right\}$ 优化。最终，综合预测模型通过递归修正实现：

${\widehat{y}}^{(L)}\left(x\right)={\widehat{y}}^{(1)}\left(x\right)+{\displaystyle \sum _{l=1}^{L-1}{\widehat{\delta }}^{(l)}}\left(x\right),$

其中${\widehat{y}}^{(1)}\left(x\right)$ 为最低精度模型预测值，修正项${\widehat{\delta }}^{(l)}\left(x\right)$ 逐层级叠加，逐步逼近高精度响应。

4. 数值实验

为验证所提出的多精度TAG过程的有效性，本节选取五种高维测试函数以及Borehole函数数值实验，并比较了四种不同方法在嵌套拉丁超立方设计生成的高低精度数据上的效果。

4.1. 高低精度测试函数

用MFTAG表示本文所提出的方法，其中变换函数为Box-Cox变换及其拓展形式。MFGP表示传统的多精度高斯过程，使用了R中的DiceKriging包[31]来实现。MFTAGH和MFGPH表示仅使用高精度数据的建模结果。其中高精度样本数量$n_h$ 和低精度样本数量$n_l$ 与维度$p$ 有关，在低维情形($p<10$ )时，$n_h=20$ ，$n_l=40$ ；当$p=10$ 时，$n_h=100$ ，$n_l=200$ ；当$p=20$ 时，$n_h=180$ ，$n_l=360$ 。每种方法重复100次，比较它们的均方预测误差(Mean Squared Prediction Error, MSE)的结果

$\frac{1}{10,000}{\displaystyle \sum _{k=1}^{10,000}{\left[{\widehat{y}}_h\left(x_i\right)-y_h\left(x_i\right)\right]}^2},$

其中$x_1,\cdots ,x_{10000}$ 由拉丁超立方设计生成，结果展示在表1中。测试函数形式如下：

Model 1来自Currinet等[28]，其中高精度函数为

$y_h=\left[1-\exp \left(-\frac{1}{2x_2}\right)\right]\frac{2300x_1^3+1900x_1^2+2092x_1+60}{100x_1^3+500x_1^2+4x_1+20}.$

低精度函数为

$\begin{array}{c}{y}_l=[y_h\left(x_1+1/20,x_2+1/20\right)+y_h\left(x_1+1/20,\max \left(0,x_2-1/20\right)\right)\\ +y_h\left(x_1-1/20,x_2+1/20\right)+y_h\left(x_1-1/20,\max \left(0,x_2-1/20\right)\right)]/4.\end{array}$

Model 2来自Cox等[32]，其中高精度函数为

$y_h=\frac{2}{3}\exp \left(x_1+x_2\right)-x_4\sin \left(x_3\right)+x_3.$

低精度函数为

$y_l=1.2y_h-1.$

Model 3来自Cox等[32]，高精度函数为

$y_h=\frac{x_1}{2}\left[\sqrt{1+\left(x_2+x_3^2\right)x_4/x_1^2}-1\right]+\left(x_1+3x_4\right)\times \exp \left[1+\sin \left(x_3\right)\right].$

低精度函数为

$y_l=\left[1+\sin \left(x_1\right)/10\right]y_h\left(x_1,x_2,x_3,x_4\right)-2x_1+x_2^2+x_3^2+\mathrm{0.5.}$

Model 4为一个加权模型，高精度函数为

$y_h={\displaystyle \sum _{j=1}^{p}j}{x}_j^2,$

低精度函数为

$y_l=1.2y_h-1.$

Model 5为Ackley函数，高精度函数为

$y_h=20\exp \left(-\frac{1}{5}\sqrt{\frac{1}{p}{\displaystyle \sum _{j=1}^p{x}_j^2}}\right)+\exp \left(\frac{1}{p}{\displaystyle \sum _{j=1}^{p}2}\pi x_j\right)-20-\exp \left(1\right),$

低精度函数为

$y_l=20\exp \left(-\frac{1}{4}\sqrt{\frac{1}{p}{\displaystyle \sum _{j=1}^p{x}_j^2}}\right)+\exp \left(\frac{1}{p}{\displaystyle \sum _{j=1}^{p}6}{x}_j\right)-20,$

Model 6为Levy函数，高精度函数为

$y_h={\sin }^2\left(\pi w_1\right)+{\displaystyle \sum _{i=1}^{d-1}{\left(w_i-1\right)}^2}\left[1+10{\sin }^2\left(\pi w_i+1\right)\right]+{\left(w_d-1\right)}^2\left[1+{\sin }^2\left(2\pi w_d\right)\right],$

低精度函数为

$y_1={\sin }^2\left(3w_1\right)+{\displaystyle \sum _{i=1}^{d-1}{\left(w_i-1\right)}^2}\left[1+10{\sin }^2\left(3w_i+1\right)\right]+{\left(w_d-1\right)}^2\left[1+{\sin }^2\left(6w_d\right)\right],$

其中$w_i=1+\frac{x_i-1}{4},\cdots ,d$ 。

Model 7为Dixon-Price函数，高精度函数为

$y_h={\left(x_1-1\right)}^2+{\displaystyle \sum _{i=2}^{d}i}{\left(2x_i^2-x_{i-1}\right)}^2$

低精度函数为

$y_l=1.15y_h-0.5\sin \left(1.5y_h\right)+0.3\frac{y_h^3}{1+{\left|y_h\right|}^2}-0.7$

Model 8为Trid函数，高精度函数为

$y_h={\displaystyle \sum _{i=1}^d }{\left(x_i-1\right)}^2-{\displaystyle \sum _{i=2}^d{x}_i{x}_{i-1}}$

低精度函数为

$y_1=1.1y_h-0.6\tanh \left(2y_h\right)+0.4{\text{e}}^{-0.5y_h^2}-\mathrm{0.8.}$

Table 1. Result of MSE

表1. MSE结果

模型	维度	MFTAG	MFGP	MFTAGH	MFGPH
Model 1	2	0.0697 (0.2539)	0.4545 (0.7038)	0.4023 (1.0663)	1.4729 (2.2093)
Model 2	4	0.0003 (0.0005)	0.0005 (0.0002)	0.0069 (0.0128)	0.0066 (0.0042)
Model 1	10	0.0062 (0.0105)	0.0459 (0.0120)	0.0258 (0.0436)	0.1431 (0.0581)
Model 2	10	0.0006 (0.0002)	0.0010 (0.0002)	0.0013 (0.0006)	0.0041 (0.0008)
Model 3	10	0.0126 (0.0056)	0.0175 (0.0027)	0.0309 (0.0094)	0.0691 (0.0150)
Model 1	20	0.0046 (0.0090)	0.0657 (0.0094)	0.0219 (0.0472)	0.1480 (0.0259)
Model 2	20	0.0006 (0.0001)	0.0038 (0.0004)	0.0016 (0.0004)	0.0089 (0.0012)
Model 3	20	0.0261 (0.0214)	0.0834 (0.0080)	0.0311 (0.0065)	0.1536 (0.0234)
Model 4	20	0.0004 (0.0001)	1.6071 (0.1631)	0.0014 (0.0002)	3.1833 (0.2911)
Model 5	20	0.4016 (0.1567)	1.9317 (0.2356)	0.4724 (0.2308)	4.5311 (0.5965)
Model 6	20	0.7832 (0.4279)	14.9139 (0.7906)	2.8917 (0.8951)	16.1598 (0.4951)
Model 7	20	0.1162 (0.0112)	0.1370 (0.0065)	0.5161 (0.1981)	0.8941 (0.1566)
Model 8	20	0.1123 (0.0574)	0.1774 (0.0198)	0.5958 (0.1566)	0.6161 (0.1561)

表1中括号内为标准差，从表中结果可以看出本文所提方法的结果在多数情形下显著优于基于传统的多精度高斯过程，表明了TAG过程在高维度的多精度计算机试验建模中的优越性。但在交互情形较为严重的情形(如Model 7和Model 8)，MFTAG有着和MFGP相近的效果，表明在该情形下TAG过程以及所使用的Box-Cox变换未能完全捕捉到交互效应，未来可采用其他变换来改进效果。

4.2. Borehole函数

接下来我们使用Borehole函数[33]来验证我们所提出的方法，即

$y_h=\frac{2\pi T_u\left(H_u-H_l\right)}{\log \left(r/r_w\right)\left[1+\frac{2L{T}_u}{\log \left(r/r_w\right)r_w^2{K}_w}+T_u/T_l\right]}.$

该函数描述了地下水通过钻孔从地表穿越两个含水层的流动过程，是一类典型的水文地质模拟问题。该函数因其非线性特性和广泛应用，是计算机试验研究中常用的测试函数之一。该模型中八个变量的取值范围分别为：$r_w=\left(0.05,0.15\right)$ ，$r=\left(100,50,000\right)$ ，$T_u=\left(63,070,115,600\right)$ ，$H_u=\left(990,1110\right)$ ， $T_l=\left(63.1,116\right)$ ，$H_l=\left(700,820\right)$ ，$L=\left(1120,1680\right)$ ，$K_w=\left(1500,15,000\right)$ 。其所对应的低精度版本来自Xiong等人[30]

$y_l=\frac{5T_u\left(H_u-H_l\right)}{\log \left(r/r_w\right)\left[1.5+\frac{2L{T}_u}{\log \left(r/r_w\right)r_w^2{K}_w}+T_u/T_l\right]}.$

在此模型中，我们选取高精度样本数量$n_h=320$ ，低精度样本数量$n_l=640$ ，每种方法仍重复100次，比较MSE结果。经计算，MFTAG、MFGP、MFTAGH、MFGPH所对应的MSE分别为0.0672 (0.0710)，0.1308 (0.0280)，0.0891 (0.0976)，0.5053 (0.0987)，其中括号内为标准差。结果表明本文所提方法仍优于传统的多精度高斯过程方法。

5. 总结

本文提出了一种基于变换加性高斯过程的多精度建模方法，用于解决高维多精度计算机试验中的建模问题。TAG过程通过变换后的响应变量建立可加结构，每个变量部分用独立的高斯过程建模。并通过引入额外的平稳高斯过程来处理可能的高阶交互效应，同时通过分步拟合和一维高斯过程初始化参数，减少了优化复杂度，适用于高维数据的建模。本文将TAG过程引入多精度计算机试验的建模中，利用嵌套设计以及层级化误差修正机制，显著提升了高维情形下的建模精度。数值实验也表明，本文所提方法在高维多精度计算机试验中效果一致优于多精度高斯过程模型，为高维多精度问题提供了一种实用解决方案。

基金项目

国家自然科学基金项目(12171462)。

NOTES

^*通讯作者。