分层Kriging模型理论相关推导

分层Kriging(Hierarchical Kriging, HK)模型是一种简单实用的变可信度代理模型，分两层或多层建立代理模型，本文以两层模型为例对分层代理模型理论进行推导。

3.1 代理模型问题的基本描述

对于一个有$m$个设计变量的优化问题，在设计空间中同时进行高可信度分析和低可信度抽样，以建立所谓变可信度模型，变可信代理模型在达到相同近似精度的条件下，可以显著提高建立代理模型的效率。设高、低可信度分析程序的抽样位置分别为：
$$\{\begin{array}{l}{\mathbf{S}}_1=[{\mathbf{x}}_1^{(1)}{\mathbf{x}}_1^{(2)}\dots {\mathbf{x}}_1^{(n_1)}{]}^T\in {\mathbb{R}}^{n_1\times m}\\ {\mathbf{S}}_2=[{\mathbf{x}}_2^{(1)}{\mathbf{x}}_2^{(2)}\dots {\mathbf{x}}_2^{(n_2)}{]}^T\in {\mathbb{R}}^{n_2\times m}\end{array}\text{\hspace{1em}(1)}$$
上式中，下标 “1” 和 “2” 分别代表高、低可信度，例如 $n_1$ 和 $n_2$ 分别代表高、低可信度样本点数（可合理假设 $n_2\gg n_1$ ）。相应的目标函数或约束函数响应值为：
$$\{\begin{array}{l}{\mathbf{y}}_1=[y_1^{(1)}{y}_1^{(2)}\dots y_1^{(n_1)}{]}^T\in {\mathbb{R}}^{n_1}\\ {\mathbf{y}}_2=[y_2^{(1)}{y}_2^{(2)}\dots y_2^{(n_2)}{]}^T\in {\mathbb{R}}^{n_2}\end{array}\text{\hspace{1em}(2)}$$
分层代理模型首先在低可信度样本数据 $({\mathbf{S}}_2,y_2)$ 的基础上建立 Kriging 模型 ${\hat{y}}_2$ ，然后以 ${\hat{y}}_2$ 为全局趋势模型，在高可信度样本数据集 $({\mathbf{S}}_1,y_1)$ 基础上建立所需的代理模型 ${\hat{y}}_1$。

3.2 HK 模型及其预估值

3.2.1 首先针对低可信度样本点建立Kriging代理模型

首先基于低可信度数据来建立低可信度代理模型以辅助预测，针对低可信度模型假设高斯静态随机过程：
$$Y_{lf}(\mathbf{x})={\beta }_{0,lf}+Z_{lf}(\mathbf{x})\text{\hspace{1em}(3)}$$
参照 Kriging 模型建立过程，可得低可信度模型表达式：
$${\hat{y}}_{lf}(\mathbf{x})={\beta }_{0,lf}+{\mathbf{r}}_{lf}^T(\mathbf{x}){\mathbf{R}}_{lf}^{-1}({\mathbf{y}}_{s,lf}-{\beta }_{0,lf}{\mathbf{F}}_1)\text{\hspace{1em}(4)}$$
上式中，${\beta }_{0,lf}=({\mathbf{F}}_1^T{\mathbf{R}}_{lf}^{-1}{\mathbf{F}}_1{)}^{-1}{\mathbf{F}}_1^T{\mathbf{R}}_{lf}^{-1}{\mathbf{y}}_{s,lf}$ ；${\mathbf{R}}_{lf}\in {\mathbb{R}}^{n_{lf}\times n_{lf}}$ 是现有的低可信度样本点组成的相关矩阵；${\mathbf{F}}_1$ 是全为1的列向量；${\mathbf{r}}_{lf}\in {\mathbb{R}}^{n_{lf}}$ 是任意点 $\mathbf{x}$ 与已知样本点间的相关矢量。

3.2.2 建立分层代理模型

不同于Cokriging 模型，将低可信度数据建立的代理模型作为趋势项，将高可信度函数视为下述高斯静态随机过程的具体实现：
$$Y(\mathbf{x})={\beta }_0{\hat{y}}_{lf}(\mathbf{x})+Z(\mathbf{x})\text{\hspace{1em}(5)}$$
将 ${\hat{y}}_{lf}(\mathbf{x})$ 乘以缩放系数 ${\beta }_0$ 作为全局趋势模型，$Z(\mathbf{x})$ 是均值为零、方差为 ${\sigma }^2$ 的静态随机过程。在设计空间不同位置处，这些随机变量存在一定的相关性（协方差），表示为：
$$Cov[Z(\mathbf{x}),Z({\mathbf{x}}^{\prime })]={\sigma }^{2}R(\mathbf{x},{\mathbf{x}}^{\prime })\text{\hspace{1em}(6)}$$
上式中， $R(\mathbf{x},{\mathbf{x}}^{\prime })$ 是相关函数，它是 $\mathbf{x}$ 和 ${\mathbf{x}}^{\prime }$ 间欧几里德距离的函数，将高可信度函数视为已知高可信度样本点的线性加权，即：
$$\hat{y}(\mathbf{x})={\mathbf{\omega }}^T{\mathbf{y}}_s\text{\hspace{1em}(7)}$$
上式中，$\mathbf{\omega }={\left[\begin{array}{ccc}{\omega }^{(1)}& {\omega }^{(2)}& \dots {\omega }^{(n)}\end{array}\right]}^T$ 是高可信度样本点的加权系数向量，用 ${\mathbf{Y}}_s=\left[\begin{array}{cccc}{Y}^{(1)}& Y^{(2)}& \dots & Y^{(n_1)}\end{array}\right]$ 代替 ${\mathbf{y}}_s=\left[\begin{array}{cccc}{y}^{(1)}& y^{(2)}& \dots & y^{(n)}\end{array}\right]$ 。寻找最优加权系数 $\mathbf{\omega }$ 使得均方差最小：
$$MSE[\hat{y}(\mathbf{x})]=E[({\mathbf{\omega }}^T{\mathbf{Y}}_s-Y(\mathbf{x}){)}^2]\text{\hspace{1em}(8)}$$
并满足无偏约束：
$$E[{\mathbf{\omega }}^T{\mathbf{Y}}_s]=E[Y(\mathbf{x})]\text{\hspace{1em}(9)}$$
对上式化简得：

以下推导过程要注意，应将 ${\hat{y}}_{lf}(\mathbf{x})$ 视为由低可信度代理模型得出的值，代入到高可信度模型后则应视为常值（此常值在设计空间不同点处的值不同）

$$\sum _{i=1}^{n_1}{\omega }^{(i)}{\hat{y}}_{lf}({\mathbf{x}}_1^{(i)})={\hat{y}}_{lf}(\mathbf{x})\text{\hspace{1em}(10)}$$

采用拉格朗日乘数法，令：
$$L(\mathbf{\omega },\mu )=E[({\mathbf{\omega }}^T{\mathbf{Y}}_s-Y(\mathbf{x}){)}^2]+\mu (\sum {\omega }^{(i)}{\hat{y}}_{lf}({\mathbf{x}}_1^{(i)})-{\hat{y}}_{lf}(\mathbf{x}))\text{\hspace{1em}(11)}$$
令：
∂ L ( ω , μ ) ∂ ω ( i ) = 2 E [ ( ω T Y s − Y ( x ) ) Y s ( x 1 ( i ) ) ] + μ y ^ l f ( x 1 ( i ) ) = 2 { ∑ j = 1 n 1 ω ( j ) [ β 0 2 y ^ l f ( x 1 ( i ) ) y ^ l f ( x 1 ( j ) ) + σ 2 R ( x 1 ( i ) , x 1 ( j ) ) ] − [ β 0 2 y ^ l f ( x ) y ^ l f ( x 1 ( i ) ) } + σ 2 R ( x 1 ( i ) , x ) ] + μ y ^ l f ( x 1 ( i ) ) = 2 σ 2 ∑ j = 1 n 1 ω ( j ) R ( x 1 ( i ) , x 1 ( j ) ) − 2 σ 2 R ( x 1 ( i ) , x ) + μ y ^ l f ( x 1 ( i ) ) = 0 (12) \begin {align} \frac {\partial L(\boldsymbol \omega,\mu)}{\partial \omega^{(i)}} &= 2E[(\boldsymbol \omega^T\mathbf Y_s-Y(\mathbf x))Y_s(\mathbf x_1^{(i)})] +\mu\hat y_{lf}(\mathbf x_1^{(i)}) \\ &= 2\{\sum_{j=1}^{n_1}\omega^{(j)}[\beta_0^2\hat y_{lf}(\mathbf x_1^{(i)})\hat y_{lf}(\mathbf x_1^{(j)})+\sigma^2R(\mathbf x_1^{(i)},\mathbf x_1^{(j)})]-[\beta_0^2\hat y_{lf}(\mathbf x)\hat y_{lf}(\mathbf x_1^{(i)})\}+\sigma^2R(\mathbf x_1^{(i)},\mathbf x)]+\mu \hat y_{lf}(\mathbf x_1^{(i)}) \\ &= 2\sigma^2\sum_{j=1}^{n_1}\omega^{(j)}R(\mathbf x_1^{(i)},\mathbf x_1^{(j)})-2\sigma^2R(\mathbf x_1^{(i)},\mathbf x)+\mu\hat y_{lf}(\mathbf x_1^{(i)}) \\ &=0 \end {align} \tag{12} ∂ω(i)∂L(ω,μ)​​=2E[(ωTYs​−Y(x))Ys​(x1(i)​)]+μy^​lf​(x1(i)​)=2{j=1∑n1​​ω(j)[β02​y^​lf​(x1(i)​)y^​lf​(x1(j)​)+σ2R(x1(i)​,x1(j)​)]−[β02​y^​lf​(x)y^​lf​(x1(i)​)}+σ2R(x1(i)​,x)]+μy^​lf​(x1(i)​)=2σ2j=1∑n1​​ω(j)R(x1(i)​,x1(j)​)−2σ2R(x1(i)​,x)+μy^​lf​(x1(i)​)=0​(12)
则：
$$\{\begin{array}{l}{\sum }_{j=1}^{n_1}{\omega }^{(j)}R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)})+\frac{\mu {\hat{y}}_{lf}({\mathbf{x}}_1^{(i)})}{2{\sigma }^2}=R({\mathbf{x}}_1^{(i)},\mathbf{x})\\ {\sum }_{i=1}^{n_1}{\omega }^{(i)}{\hat{y}}_{lf}({\mathbf{x}}_1^{(i)})={\hat{y}}_{lf}(\mathbf{x})\end{array}\text{\hspace{1em}(13)}$$
写成矩阵形式为：
$$\left[\begin{array}{cc}\mathbf{R}& \mathbf{F}\\ {\mathbf{F}}^T& 0\end{array}\right]\left[\begin{array}{c}\mathbf{\omega }\\ \hat{\mu }\end{array}\right]=\left[\begin{array}{c}\mathbf{r}\\ {\hat{y}}_{lf}(\mathbf{x})\end{array}\right]\text{\hspace{1em}(14)}$$
上式中：
$$\begin{array}{rl}& \mathbf{F}=[{\hat{y}}_{lf}({\mathbf{x}}_1^{(1)}),\dots ,{\hat{y}}_{lf}({\mathbf{x}}_1^{(n_1)}){]}^T,\hat{\mu }=\mu /(2{\sigma }^2)\\ & \mathbf{R}=(R({\mathbf{x}}_1^{(i)},{\mathbf{x}}_1^{(j)}){)}_{i,j}\in {\mathbb{R}}^{n_1\times n_1},\mathbf{r}=(R({\mathbf{x}}_1^{(i)},\mathbf{x}){)}_i\in {\mathbb{R}}^{n_1}\end{array}\text{\hspace{1em}(15)}$$
则：
$$\left[\begin{array}{cc}{\mathbf{\omega }}^T& \hat{\mu }\end{array}\right]={\left({\left[\begin{array}{cc}\mathbf{R}& \mathbf{F}\\ {\mathbf{F}}^T& 0\end{array}\right]}^{-1}\left[\begin{array}{c}\mathbf{r}\\ {\hat{y}}_{lf}(\mathbf{x})\end{array}\right]\right)}^T={\left[\begin{array}{c}\mathbf{r}\\ {\hat{y}}_{lf}(\mathbf{x})\end{array}\right]}^T{\left[\begin{array}{cc}\mathbf{R}& \mathbf{F}\\ {\mathbf{F}}^T& 0\end{array}\right]}^{-1}\text{\hspace{1em}(16)}$$
则设计空间任意未知点 $\mathbf{x}$ 处预测的响应值为：
$$\hat{y}(\mathbf{x})=\left[\begin{array}{cc}{\mathbf{\omega }}^T& \hat{\mu }\end{array}\right]\left[\begin{array}{c}{\mathbf{y}}_s\\ 0\end{array}\right]={\left[\begin{array}{c}\mathbf{r}\\ {\hat{y}}_{lf}(\mathbf{x})\end{array}\right]}^T{\left[\begin{array}{cc}\mathbf{R}& \mathbf{F}\\ {\mathbf{F}}^T& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathbf{y}}_s\\ 0\end{array}\right]\text{\hspace{1em}(17)}$$

由分块矩阵求逆公式，有：
$${\left[\begin{array}{cc}\mathbf{R}& \mathbf{F}\\ {\mathbf{F}}^T& 0\end{array}\right]}^{-1}=\left[\begin{array}{cc}{\mathbf{R}}^{-1}-{\mathbf{R}}^{-1}\mathbf{F}({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}& {\mathbf{R}}^{-1}\mathbf{F}({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}\\ ({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}& -({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}\end{array}\right]\text{\hspace{1em}(18)}$$
将上式代入式 (16) 并化简得：
$$\begin{array}{rl}\hat{y}(\mathbf{x})& =[{\mathbf{r}}^T{\mathbf{R}}^{-1}-{\mathbf{r}}^T{\mathbf{R}}^{-1}\mathbf{F}({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}+{\hat{y}}_{lf}(\mathbf{x})({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}]{\mathbf{y}}_s\\ & ={\beta }_0{\hat{y}}_{lf}(\mathbf{x})+{\mathbf{r}}^T(\mathbf{x}){\mathbf{R}}^{-1}({\mathbf{y}}_s-{\beta }_0\mathbf{F})\end{array}\text{\hspace{1em}(19)}$$
上式中： ${\beta }_0=({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}{\mathbf{y}}_s$ 是放缩系数，$V_{HK}={\mathbf{R}}^{-1}({\mathbf{y}}_s-{\beta }_0\mathbf{F})$ 只与已知样本点数据有关，可以在模型训练结束后一次性计算并存储。之后，预测任意 $\mathbf{x}$ 处的响应值只需要计算 $\mathbf{r}(\mathbf{x})$ 与 ${\mathbf{V}}_{HK}$ 间的点乘。

分层Kriging模型给出的预估值的均方差为：

M S E [ y ^ ( x ) ] = E [ ( ω T Y s − Y ( x ) ) 2 ] = E [ ( ω T Y s − Y ( x ) ) ( ω T Y s − Y ( x ) ) ] = σ 2 ( ω T R ω − 2 ω T r + 1 ) = σ 2 [ ω T ( R ω + F μ ^ ) − ω T F μ ^ − 2 ω T r + 1 ] = σ 2 [ ω T r − ω T F μ ^ − 2 ω T r + 1 ] = σ 2 [ 1 − ω T r − ω T F μ ^ ] = σ 2 [ 1 − [ ω T    μ ^ ] [ r y l f ( x ) ] ] = σ 2 [ 1 − [ r y l f ( x ) ] T [ R F F T 0 ] − 1 [ r y l f ( x ) ] ] = σ 2 { 1.0 − r T R − 1 r + [ r T R − 1 F − y ^ l f ( x ) ] ( F T R − 1 F ) − 1 [ r T R − 1 F − y ^ l f ( x ) ] T } (20) \begin {align} MSE[\hat y(\mathbf x)] & =E[(\boldsymbol \omega ^T\mathbf Y_s-\mathbf Y(\mathbf x))^2]\\ & =E[(\boldsymbol \omega ^T\mathbf Y_s-\mathbf Y(\mathbf x))(\boldsymbol \omega ^T\mathbf Y_s-\mathbf Y(\mathbf x))]\\ & = \sigma ^2(\boldsymbol \omega^T\mathbf R \boldsymbol \omega - 2\boldsymbol \omega^T\mathbf r + 1) \\ & =\sigma^2[\boldsymbol \omega^T(\mathbf R\boldsymbol \omega+\mathbf F\hat \mu)-\boldsymbol \omega^T\mathbf F\hat \mu - 2\boldsymbol \omega^T\mathbf r + 1] \\ & = \sigma ^2[\boldsymbol \omega^T\mathbf r-\boldsymbol \omega^T\mathbf F\hat\mu-2\boldsymbol \omega^T\mathbf r + 1] \\ & =\sigma ^2[1-\boldsymbol \omega^T \mathbf r-\boldsymbol \omega^T \mathbf F\hat\mu] \\ &=\sigma ^2[1-[\boldsymbol \omega ^T \ \ \hat \mu]\begin {bmatrix}\mathbf r \\ \boldsymbol y_{lf}(\mathbf x) \end {bmatrix}]\\&= \sigma ^2[1-\begin {bmatrix} \mathbf r \\ y_{lf}(\mathbf x)\end {bmatrix}^T\begin{bmatrix} \mathbf R &\mathbf F \\ \mathbf F^T & 0 \end {bmatrix}^{-1}\begin {bmatrix}\mathbf r \\ \boldsymbol y_{lf}(\mathbf x) \end {bmatrix}] \\ &=\sigma^2 \{ 1.0-\mathbf r^T\mathbf R^{-1}\mathbf r+[\mathbf r^T\mathbf R^{-1}\mathbf F-\hat y_{lf}(\mathbf x)](\mathbf F^T\mathbf R^{-1}\mathbf F)^{-1}[\mathbf r^T\mathbf R^{-1}\mathbf F-\hat y_{lf}(\mathbf x)]^T \} \end {align} \tag {20} MSE[y^​(x)]​=E[(ωTYs​−Y(x))2]=E[(ωTYs​−Y(x))(ωTYs​−Y(x))]=σ2(ωTRω−2ωTr+1)=σ2[ωT(Rω+Fμ^​)−ωTFμ^​−2ωTr+1]=σ2[ωTr−ωTFμ^​−2ωTr+1]=σ2[1−ωTr−ωTFμ^​]=σ2[1−[ωT  μ^​][rylf​(x)​]]=σ2[1−[rylf​(x)​]T[RFT​F0​]−1[rylf​(x)​]]=σ2{1.0−rTR−1r+[rTR−1F−y^​lf​(x)](FTR−1F)−1[rTR−1F−y^​lf​(x)]T}​(20)

分层Kriging 模型的相关函数，模型训练过程与Kriging模型十分相似，可以参考之前的博客文章(链接: link)，在此不再赘述。

参考文献

[1] HAN Z H , Görtz S. Hierarchical Kriging Model for Variable-Fidelity Surrogate Modeling[J]. AIAA Journal, 2012, 50(9):1885-1896.