多精度代理模型简明教程

引言

本文将简要介绍kriging代理模型和多层代理模型的原理及使用方法，以及对其程序进行简要说明。

目前随着计算能力的不断提升，我们已经能够对一些物理过程进行精细建模和高精度计算，但是随着计算规模的提升，高精度计算带来的时间和计算成本也常常难以让人忍受。而使用代理模型可以有效减少计算量，同时获取我们感兴趣的参数。代理模型亦称为响应面模型或元模型，使用代理模型进行设计优化已经成为重要的一种优化方式。

与其他简单形式的代理模型，例如多项式模型和径向基函数模型等相比，kriging模型具有如下两个特点：1）预测精度高，2）可以给出预测准确度估计。因此，在实际应用中kriging模型正变得越来越受重视。许多文章和书籍都将kriging模型描述为一种“modelling the function as a realization of a stochastic process”，这是一种基于最小化期望平方误差（minimize the expected squared prediction error）的模型，并且满足预测结果是无偏估计，而且是其他已知结果$y_i$的线性组合。通过一系列初始数据可以得到kriging代理模型的参数，从而能够构建代理模型，这一过程称为训练，而这些初始数据称为训练集。

代理模型的预测精度取决于其训练集的精度，可以想象训练集的精度越高、数量越多，训练得到的代理模型精度也会越高。但是在实际应用过程中，如果要获取大量高精度训练集需要耗费过多的时间和计算资源，因此为了解决这一问题，我们可以考虑使用多精度代理模型（multi-fidelity surrogate model）。多精度代理模型的出发点是希望能够用大量低精度的数据和少量高精度的数据构建高精度代理模型，这样就可以大大江都获取高精度代理模型的难度。

接下来我们将简要介绍kriging代理模型、多精度kriging代理模型的原理和推导过程，详细内容请参见如下两篇文献：

• kriging代理模型：A Taxonomy of Global Optimization Methods Based on Response Surfaces by Donald R. Jones.
• 多精度kriging代理模型：Multi-fidelity optimization via surrogate modelling by Alexander I. J. Forrester, Andras Sobester and Andy J. Keane

1. kriging代理模型简介

1.1 概述

在工程界我们常常会遇到这样一些问题，希望通过建立输入参数-输出参数之间的映射关系，使得给定新的输入来得到新的输出，即预测值。根据输出值的类型为连续型和离散型可以分为回归问题和分类问题，处理回归问题的一种重要方法是高斯过程回归(Gaussian Process Regression)，kriging代理模型是高斯过程回归方法的一种应用。回归问题可以有如下数学描述：

假设有训练集$ D= {(\boldsymbol{x}_i,y_i)|i = 1,2,\cdots,n} = (\boldsymbol{X},\boldsymbol{y}). $其中：$\boldsymbol{x}_i \in R^d$为d维输入矢量，$\boldsymbol{X} = [\boldsymbol{x}_1, \boldsymbol{x}2, \cdots, \boldsymbol{x}n]$ 为$d\times n$维矩阵，$y_i \in R$为相应的输出标量，$\boldsymbol{y}$为输出矢量。回归的任务是希望通过训练学习建立输入$\boldsymbol{X}$与输出$\boldsymbol{y}$之间的之间的映射关系$(f(\cdot) : R^d \mapsto R)$，预测出新测试点$\boldsymbol{x}*$对应的最有可能的输出值 $y*.$

高斯过程是任意有限个随机变量具有联合高斯分布的集合，其性质完全由均值函数$m(\boldsymbol{x})$和协方差函数$K(\boldsymbol{x},\boldsymbol{x'})$决定，定义如下：
$$
m(\boldsymbol{x}) = E[f(\boldsymbol{x})] \
K(\boldsymbol{x},\boldsymbol{x'}) = E[(f(\boldsymbol{x}) - m(\boldsymbol{x})(f(\boldsymbol{x'}) - m(\boldsymbol{x'})]
$$
接下来我们简要推导一下kriging模型的实现方法，以对该方法有一些感性的理解。kriging模型为：
$$
Z(s) = \mu + \epsilon(s)
$$
其中，等式左侧为观测值，等式右侧分别为平均值和误差，在某一空间维度中表示如下1：

image

上面这种模型被称为普通kriging模型，因为我们在公式中引入了常数平均值。在很多情况下，平均值函数还可以表示为一次函数、二次函数等不同种类，但是简单起见我们这里只讨论平均值为常数的情况，下文中的kriging模型均指普通kriging模型。

1.2 具体方程推导

1.2.1 求解超参数

接下来我们推导一下kriging模型的具体表达形式。假设我们已知一系列观测点${(\boldsymbol{x}_i,y_i)|i = 1,2,\cdots,n} $，同时我们对这些观测点的取值具有一定的不确定度。为了度量这些不确定度，我们假设单个观测点的分布服从高斯分布，即$y\sim N(\mu,\sigma^2)$，即我们假设这些观测值有很大概率落在区间$[\mu-3\sigma,\mu+3\sigma]$内。现在考虑这些观测点的协方差函数为$K(X,X)$，这是一个$n\times n$的矩阵。这样我们可以得到这$n$个观测点的联合高斯分布为：
$$
\boldsymbol{y} \sim N(\mu, K(X,X))
$$
现在我们考虑协方差函数$K(X,X)$的具体表达形式。由统计知识，
$$
K(X,X) = \sigma^2\boldsymbol{R}
$$
$\boldsymbol{R}$为$n\times n $的相关系数矩阵，里面的元素为不同观测点之间的相关系数，接下来我们求相关系数矩阵的具体表达形式。现在考虑两个观测点$\boldsymbol{x}_i$和$\boldsymbol{x}j$，如果我们的预测函数是连续的，那么当这两个观测点之间的距离$||\boldsymbol{x}i-\boldsymbol{x}j||$之间的距离很小，这两个观测点的观测值$y(\boldsymbol{x_i})$和$y(\boldsymbol{x_j})$也会很接近，即这两个观测值对应的随机变量之间高度相关。特别的，我们可以用高斯形式的相关系数方程来描述这样的关系：
$$
Corr[y(\boldsymbol{x_i}),y(\boldsymbol{x_j})] = \text{exp}\bigg(-\sum{l=1}^d \theta_l|\boldsymbol{x}{il}-\boldsymbol{x}{jl}|^{{pl}}\bigg)
$$
其中，$d$为观测点$\boldsymbol{x}$的维数。如果两个观测点之间的距离很近，那么它们的相关系数就会趋近于1，反之则趋近于0.两个参数$\theta_l$和$p_l$比较关键，其中$\theta_l$可以理解为相关系数对第$l$个维度的坐标距离的敏感程度，$\theta_l$越大，相关系数就对该维度的距离越敏感，下降得也越快，反之则下降得越慢。$p_l$可以衡量预测函数的光滑程度，$p_l$越接近$2$，预测函数就越光滑，可以处理一些梯度问题，如果越接近$1$，预测函数就越粗糙，不利于处理梯度问题。

经过以上分析，我们可以通过$\mu,\sigma^2,\theta_l,p_l(l=1,2,\cdots,d)$这几个参数确定这些观测点的联合高斯分布。所以现在需要求解这些参数来进行之后的预测工作，我们的思路是通过求解最大似然函数来求解参数值。当我们已知这些观测点的观测值时，意味着我们可以求出他们的似然函数表示如下：
$$
\frac{1}{(2\pi)^{n/2}(\sigma2)^{{n/2}|\boldsymbol{R}|}{1/2}}\text{exp}\bigg[\frac{-(y-\boldsymbol{1}\mu)'\boldsymbol{R}^{{-1}(y-\boldsymbol{1}\mu)}{2\sigma}2}\bigg]
$$
通过选取这些参数使得似然函数最大化意味着我们得到的这些观测点和我们的预测模型符合得最好。接下来我们需要对上式取对数，然后分别对$\mu$和$\sigma^2$进行求导，可以得到这两个参数的最佳观测值，它们都是关于相关系数矩阵$\boldsymbol{R}$的函数：
$$
-\frac{n}{2}\text{log}(\widehat{\sigma}^{2)-\frac{1}{2}\text{log}(|\boldsymbol{R}|)-\frac{(\boldsymbol{y-1}\mu)'\boldsymbol{R}}{-1}(\boldsymbol{y-1}\mu)}{2\sigma^2}\
\widehat{\mu} = \frac{\boldsymbol{1'R^{{-1}y}}{\boldsymbol{1'R}{-1}1}}\
\widehat{\sigma}^2 = \frac{(\boldsymbol{y-1}\widehat{\mu})'\boldsymbol{R}^{-1}(\boldsymbol{y-1}\widehat{\mu})}{n}
$$
将这两个参数代入$\text{log}$形式的似然函数，得到所谓“concentrated log-likelihood function"，忽略常数项得到后得到：
$$
-\frac{n}{2}\text{log}(\widehat{\sigma}^2)-\frac{1}{2}\text{log}(|\boldsymbol{R}|)
$$
该方程仅依赖与相关系数矩阵$\boldsymbol{R}$，而相关系数矩阵依赖于参数$\theta_l$和$p_l$。这样，我们的求解思路就很明确了。首先通过最优化方法寻找使concentrated log-likelihood function 最大化的参数$\theta_l$和$p_l$，可以使用遗传算法或粒子群算法等随机类优化算法避免陷入局部最优。然后，代入求解参数$\mu$和$\sigma^{2$的估计值$\widehat{\mu}$和$\widehat{\sigma}}2$。

1.2.2 预测值与误差估计

到此为止我们可以得到我们想要的这些参数，接下来我们分析一下预测函数的实现方式。简单来说，假设我们现在有一个新的测量点$\boldsymbol{x}^*$，$y$是该测量点的预测值，那么根据这$(n+1)$个测量点我们就得到了一个新的似然方程，我们的目标就是要寻找$y^$，使得新的似然方程最大化。求得的$y^*$也就是我们希望得到的kriging预测点。

假设$\widetilde{\boldsymbol{y}}=(\boldsymbol{y}'\quad y^{*)$为加入新预测值$y}$的增广预测值向量，$\boldsymbol{r}$代表新预测值与已知预测值的相关系数向量。这样我们得到新的$(n+1)\times(n+1)$增广协方差矩阵$\widetilde{R}$:
$$
\widetilde{R} = \left( \begin{array}{cc}
\boldsymbol{R}&\boldsymbol{r}\
\boldsymbol{r}'&1\
\end{array}
\right)
$$
回忆之前得到的$\text{log}$形式的似然函数，如果用增广协方差矩阵代替原来的协方差矩阵，那么该似然函数依赖于$y^$的部分为：
$$
-\frac{(\widetilde{\boldsymbol{y}}-\boldsymbol{1}\widehat{\mu})'\widetilde{\boldsymbol{R}}^{{-1}(\widetilde{\boldsymbol{y}}-\boldsymbol{1}\widehat{\mu})}{2\widehat{\sigma}}2}
$$
这一部分主要难点是如何求解增广协方差矩阵的逆矩阵$\widetilde{\boldsymbol{R}}^{{-1}$，幸运的是已经有人帮我们求好了这一部分，由于形式较复杂，这里就不展出了。由于我们希望求解$y}$，使得增广似然方程中依赖于$y^$的部分最大化，所以我们将上式具体展开，并对$y^$求导，令导数等于零，得到如下方程：
$$
\bigg\frac{-1}{\widehat{\sigma}^{2(1-\boldsymbol{r}'\boldsymbol{R}}{-1}\boldsymbol{r})}\bigg+\bigg[\frac{\boldsymbol{r'\boldsymbol{R}^{{-1}}(\boldsymbol{y-1}\widehat{\mu})}{\widehat{\sigma}}2(1-\boldsymbol{r}'\boldsymbol{R}^{-1}\boldsymbol{r})}\bigg]=0
$$
由上式可以求得kriging预测值$\widehat{y}(\boldsymbol{x^})$:
$$
\widehat{y}(\boldsymbol{x^{*})=\widehat{\mu}+r'\boldsymbol{R}}{-1}(\boldsymbol{y-1}\widehat{\mu})
$$
进一步地，我们可以求出在测量点$\boldsymbol{x}^$处的误差为：
$$
s^{2(\boldsymbol{x}}) = \widehat{\sigma}^{2\bigg[1-\boldsymbol{r}'\boldsymbol{R}}{-1}\boldsymbol{r} + \frac{(1-\boldsymbol{r}'\boldsymbol{R}^{{-1}\boldsymbol{r})}2}{\boldsymbol{1}'\boldsymbol{R}^{-1}\boldsymbol{1}}\bigg]
$$
该误差方程等号右端最右侧一项可以理解为对$\mu$预测的不确定性引起的预测值的不确定性。值得注意的是，在已知测量点上，该误差为零，这一结论可以根据误差函数严格证明，这里我们就不再详述。

在我们进行数据输入时，为了消除单位变动对距离的影响，我们常常将输入数据先进行归一化，然后在计算结束的时候再复原数据。

2. 多精度代理模型cokriging简介

2.1 求解超参数

之前已经在序言里阐述过，多精度代理模型的主要思想是希望能够用计算量较低的低精度数据代替大部分需要消耗大量计算资源的高精度数据，并达到和高精度代理模型接近的计算精度，从而减少高精度数据的计算量。假设我们有两套数据，分别是在观测点$\boldsymbol{X_e}$的高精度观测值$\boldsymbol{y}_e$和在观测点$\boldsymbol{X}_c$的低精度观测值$\boldsymbol{y}_c$，这里最好满足$(\boldsymbol{X}_e\subset\boldsymbol{X}_c)$.我们将这些观测点组合起来得到联合观测点：
$$
\boldsymbol{X} = \left(\begin{array}{c}
\boldsymbol{X}_c\
\boldsymbol{X}_e\
\end{array}\right)
=(\boldsymbol{x}_c^{1,\cdots,\boldsymbol{x}_c}{n_c},\boldsymbol{x}_e^{1,\cdots,\boldsymbol{x}_e}{n_e})^T
$$
在kriging模型中，在观测点$\boldsymbol{X}$处的预测值是一个高斯随机变量，因此这里我们有联合预测值
$$
\boldsymbol{Y} = \left(\begin{array}{c}
\boldsymbol{Y}_c(\boldsymbol{X}_c)\
\boldsymbol{Y}_e(\boldsymbol{X}_e)\
\end{array}\right)
=(Y_c(\boldsymbol{x}_c^{1),\cdots,Y_c(\boldsymbol{x}_c}{n_c}),Y_c(\boldsymbol{x}_e^{1),\cdots,Y_c(\boldsymbol{x}_e}{n_e}))^T
$$
这里我们假设在观测点$\boldsymbol{x}_i$处的预测值具有Markov性质，即$cov{Y_e(\boldsymbol{x}^{i),Y_c(\boldsymbol{x})|Y_c(\boldsymbol{x}}i)}=0,\forall\boldsymbol{x}\ne\boldsymbol{x}^i$^[1]，这个性质在我们求解cokriging的协方差矩阵时有很大帮助。高斯过程$Z_c(\cdot)$和$Z_e(\cdot)$代表低精度和高精度kriging模型的当地特征，它们之间的关系表示如下：
$$
Z_e(\boldsymbol{x}) = \rho Z_c(\boldsymbol{x}) + Z_d(\boldsymbol{x}).
$$
其中，$\rho$是一个标量，$Z_d(\boldsymbol{x})$代表高斯过程$Z_e(\boldsymbol{x})$和$\rho Z_c(\boldsymbol{x})$之间的差别。这样，我们就可以求出cokriging模型的协方差矩阵$cov{Y(\boldsymbol{x}),Y(\boldsymbol{x})} $：
$$
\begin{array}{l}
cov{\boldsymbol{Y}_c(\boldsymbol{X}_c),\boldsymbol{Y}_c(\boldsymbol{X}_c)} &=
cov{Z_c(\boldsymbol{X}_c),Z_c(\boldsymbol{X}_c)}\
&= \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_c,\boldsymbol{X}_c)\

cov{\boldsymbol{Y}_e(\boldsymbol{X}_e),\boldsymbol{Y}_c(\boldsymbol{X}_c)} &=
cov{\rho Z_c(\boldsymbol{X}_e) + Z_d(\boldsymbol{X}_e),Z_c(\boldsymbol{X}_c)}\
&=\rho \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_c)\

cov{\boldsymbol{Y}_e(\boldsymbol{X}_e),\boldsymbol{Y}_e(\boldsymbol{X}_e)} &=
cov{\rho Z_c(\boldsymbol{X}_e) + Z_d(\boldsymbol{X}_e),\rho Z_c(\boldsymbol{X}_e) + Z_d(\boldsymbol{X}_e)}\
&= \rho^2\sigma_c2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_e) + \sigma_d^2\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)
\end{array}
$$
这样，我们就可以得到更为复杂的协方差矩阵：
$$
\boldsymbol{C} = \left(\begin{array}{cc}
\sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_c,\boldsymbol{X}_c) &
\rho \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_c,\boldsymbol{X}_e)\
\rho \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_c)&
\rho^2\sigma_c2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_e) + \sigma_d^2\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)\
\end{array}\right)
$$
从协方差矩阵的构成元素中，我们可以清楚地看到这里我们需要求解更多的超参数，分别是：$\boldsymbol{\theta}_c,\boldsymbol{\theta}_d,\boldsymbol{p}_c,\boldsymbol{p}_d$和标量$\rho$。由上一章节的kriging模型，我们可以根据低精度数据求出$\boldsymbol{\theta}_c$和$\boldsymbol{p}_c$，接下来我们求解剩下的参数。首先我们由高斯过程$Z_c(\cdot)$和$Z_e(\cdot)$之间的关系，定义一个新的参数：
$$
\boldsymbol{d} = \boldsymbol{y}_e - \rho\boldsymbol{y}_c(\boldsymbol{X}_e)
$$
其中$\boldsymbol{y}_c(\boldsymbol{X}_e)$在观测点$\boldsymbol{x}_e$处的低精度观测值，一般情况下我们选择的高精度观测点是低精度观测点的子集，所以能够直接得到该值。如果某些高精度观测点不属于低精度观测点集，那么我们用构造的低精度kriging模型预测值来求出$\boldsymbol{y}_c(\boldsymbol{X}_e)$。为了求出$\boldsymbol{\theta}_d,\boldsymbol{p}_d,\rho$，我们构造出参数$\boldsymbol{d}$的$\text{log}$形式的似然函数：
$$
-\frac{n_e}{2}\text{log}(\widehat{\sigma}_d^{2)-\frac{1}{2}\text{log}(|\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)|)-\frac{(\boldsymbol{d-1}\mu_d)'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)}{-1}(\boldsymbol{d-1}\mu_d)}{2\sigma_d^2}
$$
类似地，通过求解最大似然函数，分别对$\mu_d$和$\sigma^2_d$进行求导，可以求出：
$$
\begin{array}{l}
\widehat{\mu}_d &= \frac{\boldsymbol{1'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)^{{-1}d}}{\boldsymbol{1'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)}{-1}1}}\
\widehat{\sigma}^2_d &= \frac{(\boldsymbol{d-1}\widehat{\mu}_d)'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)^{-1}(\boldsymbol{d-1}\widehat{\mu}_d)}{n_e}
\end{array}
$$
将上式代入似然函数，可以得到新的似然函数：
$$
-\frac{n_e}{2}\text{log}(\widehat{\sigma}_d^2)-\frac{1}{2}\text{log}(|\text{det}(\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e))|)
$$
通过最大化求解上面的似然函数，就可以得到我们需要的参数$\boldsymbol{\theta}_d,\boldsymbol{p}_d,\rho$，同样我们可以使用遗传算法或者粒子群算法等随机类算法来进行求解，这样可以避免陷入局部最优点。当然如果观测点的数量和维度很大，为了减少计算量我们也可以假设参数$\boldsymbol{\theta}_c,\boldsymbol{\theta}_d$均为标量，但是这样可能会降低模型的精确度。

2.2 预测值与误差估计

求出超参数之后，我们可以给出cokriging模型的预测方程：
$$
\widehat{y}_e(\boldsymbol{x}^{n_e+1}) = \widehat{\mu} + \boldsymbol{c}^{T\boldsymbol{C}}{-1}(\boldsymbol{y-1}\widehat{\mu})
$$
其中，
$$
\boldsymbol{c} = \left(\begin{array}{c}
\widehat{\rho}\widehat{\sigma}_c^{2R_c(\boldsymbol{X}_c,\boldsymbol{x}}{n+1})\
\widehat{\rho}\widehat{\sigma}_c^{2R_c(\boldsymbol{X}_e,\boldsymbol{x}}{n+1}) + \widehat{\sigma}_d^{2R_d(\boldsymbol{X}_e,\boldsymbol{x}}{n+1})
\end{array}\right)
$$
该向量是新的预测点和已知的低精度和高精度观测点的相关系数矩阵。平均值常数$ \widehat{\mu} = \frac{\boldsymbol{1'C^{{-1}y}}{\boldsymbol{1'C}{-1}1}}$ 。

进一步地，我们还可以给出在预测点的误差估计：
$$
s^2(\boldsymbol{x}) = \widehat{\rho}^{2\widehat{\sigma}_c}2 + \widehat{\sigma}_d^2 - \boldsymbol{c^TC{-1}c}
$$
如果新的预测点$\boldsymbol{x}^{n_e+1}\in\boldsymbol{X}_e$，那么该点的误差估计为$0$，可以根据上式严格证明，这里我们就不再详述。

3. 多精度代理模型示例

接下来，我们希望通过一个例子来阐述一下cokriging模型的具体实现效果。假设我们现在分别有高精度模型$f_e$和低精度模型$f_c$，表达式为：
$$
f_e(x) =(6x-2)^2\sin(12x-4),x\in[0,1]\
f_c(x) = 0.5(6x-2)^2\sin(12x-4)+10(x-0.5)+5,x\in[0,1]
$$
分别取低精度模型和高精度模型的观测点$\boldsymbol{X}_c = {0.0,0.1,0.2,0.3,0.4，0.5，0.6，0.7，0.8，0.9，1.0}$和$\boldsymbol{X}_e = {0.0,0.4,0.6,1.0}$，值得注意的是高精度模型的观测点包含低精度模型观测点的两个端点。利用这些观测点和观测值训练cokriging模型，然后取得剩余区间内的高精度预测值曲线，并与真实的函数曲线进行对比，考察cokriging模型的准确度。结果如下图所示：

图片中绿色的点为低精度观测值。可以看到，cokriging的预测曲线与真实函数曲线几乎重合，这说明我们的cokriging模型效果很好，仅仅用了四个高精度观测值就取得了很好的预测结果。

误差曲线如下图所示。可以看到误差在四个高精度观测点为零，最大方差在0.001左右。

这里也给出模型参数以供参考：
$$
\theta_c = 15.715,p_c =1 .9998\
\theta_d = 0.1,p_d = 1.99985,\rho=1.9836
$$

---

1. Kennedy M C, O'Hagan A. Predicting the output from a complex computer code when fast approximations are available[J]. Biometrika, 2000, 87(1):1-13. ↩