(十三)Kriging回归

https://zhuanlan.zhihu.com/p/377620800

一、前言

克里金(Kriging)模型是贝叶斯优化的基础，贝叶斯优化在如今的工程中应用得非常广泛。我自己的研究方向也跟克里金模型有关，最近一直在研究克里金模型是如何推导出来的。

在看文献的过程中，我非常疑惑。因为从文献来看，克里金模型的推导似乎有两种方式，但是两种方式推导出来的结果好像又不完全相同。所以我每次在使用Kriging模型的时候老是觉得不太放心。今天突然看到一篇文章，有详细的推导过程。我又把两种方式仔细看了一遍，终于发现了两者的共同点。另外一点是网上对于高斯回归的推导很多，对于克里金模型的推导却很少，让很多人不太能理解为什么高斯回归是克里金是一个东西，两者从形式上看相似，推导过程却完全不同。所以，今天在这里跟大家一起分享一下。

二、克里金模型的两种形式

形式一[1]：

文献里最常出现的形式：

{y^1(x)=μ^+rTC−1(y−1μ^)s12(x)=σ^2[1−rTC−1r+(1−1TC−1r)21TC−11]

{μ^=1TC−1y1TC−11σ^2=(y−1μ^)TC−1(y−1μ^)n

我们假设有n个数据点，1个预测点。其中 y^(x) 和s2(x)就是我们所需要的预测点x的均值和方差。 r是数据点X和预测点x之间的协方差矩阵。C是数据点X之间的协方差矩阵。y是数据点的目标值。1是n*1的矩阵。 和μ^(x)和σ^(x) 是为了求得最终结果的中间变量。

形式二[2]：

DACE是MATLAB里面实现Kriging的工具包，里面的PDF说明了Kriging模型是如何推导过来的。而且提到Kriging的几乎所有文献在具体推导时都会让读者去看1989年的这篇文献。

{y^2(x)=f(x)Tβ∗+rT(x)γ∗s22(x)=σ2(1+uT(FTR−1F)−1u−rTR−1r)

其中其中{β∗=(FTR−1F)−1FTR−1YRγ∗=Y−Fβ∗u=FTR−1r−f(x)σ2=1m(Y−Fβ∗)TR−1(Y−Fβ∗)

里面的R与形式一中的C定义相同，是数据点之间的协方差。Y是数据点的目标值。F是数据点的基函数矩阵，f(x)是未知点的基函数矩阵。

[1]Zhan, D., et al. (2017). "Expected Improvement Matrix-Based Infill Criteria for Expensive Multiobjective Optimization." IEEE Transactions on Evolutionary Computation 21(6): 956-975.

[2]Sacks, Jerome, et al. "Design and analysis of computer experiments."Statistical science(1989): 409-423.

三、二者的联系

从形式上来看，两种方式的定义完全不同，看不出来有什么联系。实际上两种形式是有内在联系的。

形式一是形式二在基函数为0次函数(即f(x)=1)时的特例

我们下面给出证明：

我们只考虑预测点为一个点的情况。当基函数是0次函数,f(x)=1时，F=1是n*1的矩阵，f=1是一个标量。

β∗=(1TR−11)−11TR−1Y=1TR−1Y1TR−11=μ^

γ∗=R−1(Y−1μ^)

将上面两式带入 y^2(x) ，得

y^2(x)=1β∗+rTγ∗=μ^+rTR−1(Y−1μ^)=y^1(x)

s22(x)=σ2(1+uTu1TR−11−r−1RTr)=σ2(1+‖u‖21TR−11−r−1RTr)=σ2(1+(1−1TR−1r)21TR−11−r−1RTr)=s12(x)

综上所述：两种形式的Kriging模型是等价的。形式一是形式二的特例。

四、推导过程

下面我们详细介绍两种形式的推导过程。形式一我们会推导至预测值，方差(MSE)的推导可以从形式二中得出。形式二我们会详细推导如何计算预测值和方差。

形式一：

{y^1(x)=μ^+rTC−1(y−1μ^)s12(x)=σ^[1−rTC−1r+(1−1TC−1r)21TC−11]

{μ^=1TC−1y1TC−11σ^2=(y−1μ^)TC−1(y−1μ^)n

形式二:

{y^2(x)=f(x)Tβ∗+rT(x)γ∗s22(x)=σ2(1+uT(FTR−1F)−1u−r−1Rr)

其中其中{β∗=(FTR−1F)−1FTR−1YRγ∗=Y−Fβ∗u=FTR−1r−f(x)σ2=1m(Y−Fβ∗)TR−1(Y−Fβ∗)

其中的f(x)是回归的模型，可以选择0次，1次，2次多项式，r(x)是相关函数，可以选择高斯核函数，指数核函数等等。

4.1 形式一模型建立

证明：

对于给定的数据集 X={x1,x2,…,xn}T ,所对应的目标函数为 y={y1,y2,…,yn}T 。克里金法假设所有数据之间都服从n维的正态分布。所以目标函数 y 是一个随机过程，里面的每一个变量 yi 都是一个随机变量。我们把这个随机过程记做：

(Y(x1)Y(x1)⋮Y(xn))∼N(μ,C)

多维随机过程里面最重要的是它的均值和协方差。我们把均值取为常数 1μ,1 是n*1的矩阵。协方差取为

cor[Y(xi),Y(xl)]=exp(−∑j=1kθj|xji−xjl|2)

C=(cor(Y(x1),Y(x1)),…,cor(Y(x1),Y(xn))⋮,⋱,⋮cor(Y(xn),Y(x1)),…,cor(Y(xn),Y(xn)))

定义了均值和方差，我们就可以写出，Y在我们定义的 和μ和σ 下，某一特定的值下的条件概率为

L(Y1,Y2,…,Yn|μ,σ)=1(2πσ2)n/2exp(−Σ(Yi−μ)22σ2)

这个条件概率可以表示为

L=1(2πσ2)n/2|C|1/2exp[−(y−1μ)TC−1(y−1μ)2σ2]

为了简化，我们对L取一个对数

ln(L)=−n2ln(2π)−n2ln(σ2)−12ln|C|−(y−1μ)TC−1(y−1μ)2σ2

有了上面的式子，问题就来了，我们应该如何选取 和μ和σ 使得Y出现在某一特定值的概率达到最大呢？答案就是对上式分别求 和μ和σ 的偏导，使得导函数为0。

我们先对 μ 求导，此时其他参数都是常数，所以我们只需要关心ln(L)的最后一项。

∂lnL∂μ=−12σ2∂(y−1μ)TC−1(y−1μ)∂μ=0⟺∂(y−1μ)TC−1(y−1μ)∂μ=0

下面问题就变成了对上面的二次矩阵求偏导，为了更加清晰地展示求导过程，我们不妨设

C−1={cij,i,j=1,…,n}

所以

(y−1μ)TC−1(y−1μ)=∑i=1n∑j=1n(yi−μ)(yj−μ)cij=∑i=1n∑j=1n(μ2−μ(yi+yj)+yiyj)cij

∂(y−1μ)TC−1(y−1μ)∂μ=∑i=1n∑j=1n(2μ−(yi+yj))cij=0

∴2μ∑i=1n∑j=1ncij=∑i=1n∑j=1ncijyi+∑i=1n∑j=1ncijyj

因为C是协方差矩阵，对称矩阵，所以C的逆矩阵也是对称矩阵。

∴2μ∑i=1n∑j=1ncij=∑j=1n∑i=1n1jcjiyi+∑i=1n∑j=1n1icijyj⟺2μ(1TC−11)=1TC−1y+1TC−1y

所以使得出现概率最大的 μ 是

μ^=1TC−1y1TC−11

我们再对 σ 求导，把 μ^ 带入得：

∂lnL∂σ=−n22σσ2−−22σ3(y−1μ^)TC−1(y−1μ^)=0

∴σ^2=(y−1μ^)TC−1(y−1μ^)n

把 和μ^和σ^ 带入lnL，忽略掉常数项，我们可以得到

lnL≈−n2ln(σ^2)−12ln(|C|)

这个式子里面的C有超参数 θ ，因此超参数的调节就是选取合适的 θ 使得lnL最大，可以用遗传算法或者是多初始点算法求得。

4.2 形式一预测

前面的工作中我们利用最大似然概率的方法得到了先验参数 和μ^,σ^和θi 的值，下面我们就详细推导如何利用这些参数来对未知点进行预测。

我们把观察到的数据记做 y ,需要预测的点(这里只考虑一个点)为 y^ ，我们将观察到的数据和要预测的点放到一起 y~={yT,y^}T ,对应的协方差矩阵也发生了改变：

C~=(CrrT1)

对应的对数似然概率为

ln(L)=−n2ln(2π)−n2ln(σ^2)−12ln|C~|−(y~−1μ^)TC~−1(y~−1μ^)2σ^2

现在的问题实际上变成了 y~ 到底取多少能够让ln(L)达到最大。因为ln(L)中前三项都与 y~ 无关，我们只需要考虑最后一项。

lnL≈−(y−1μ^y^−μ^)T(CrrT1)−1(y−1μ^y^−μ^)2σ^2

求逆是一件很麻烦的事情，下面要做的是如何把中间的逆矩阵表示出来，利用部分求逆的方法，我们可以得到

C~−1=(ABDE)=(C−1+C−1r(1−rTC−1r)−1rTC−1−C−1r(1−rTC−1r)−1−(1−rTC−1r)−1rTC−1(1−rTC−1r)−1)

详细过程为：

lnL=−(y−1μ^)TA(y−1μ^)+(y−1μ^)TB(y^−μ^)+(y^−μ^)TD(y−1μ^)+(y^−μ^)TE(y^−μ^)2σ^2

把与 y^ 无关的项都去掉，利用 B=DT ，得：

lnL≈−−12σ^2(1−rTC−1r)(y^−μ^)2+rTC−1(y−1μ^)σ^2(1−rTC−1r)(y^−μ^)

∂lnL∂y^=−−1σ^2(1−rTC−1r)(y^−μ^)+rTC−1(y−1μ^)σ^2(1−rTC−1r)=0

所以

y^1(x)=μ^+rTC−1(y−1μ^)

得证！

4.3 形式二的推导

4.3.1 问题定义

已知数据输入为 S=[s1,s2,..,sm]T , si∈Rn，输出数据为 Y=[y1,y2,…,ym]T , yi∈R1。矩阵表示为：

S=(s11,s12,…,s1ns21,s22,…,s2n⋮,⋮,⋱,⋮sm1,sm2,…,smn) , Y=(y1y2⋮ym)

输入和输出都满足(0,1)的高斯分布：

μ[S:j]=0,Var[S:iS:j]=1,j=1,2,…,n

，μ(Y)=0，Var[YY]=1

考虑一个线性模型：

y^(x)=∑i=1pβifi(x)+z(x)=[f1(x),f2(x),…,fp(x)]β+z(x)=f(x)Tβ+z(x)

其中,p代表基函数的个数,β是基函数的权重。z(x)是一个随机过程。

E(z(w)z(x))=σ2R(θ,w,x)

4.3.2 预测的推导

Kriging的总体想法是希望用已知点函数值的加权求和来表示位置点的函数值。我们先给出如下矩阵的定义，方便推导：

已知点的基函数矩阵 F=[f(s1),…,f(sm)]T

F=(f1(s1),f2(s1),…,fp(s1)f1(s2),f2(s2),…,fp(s2)⋮,⋮,⋱,⋮f1(sm),f2(sm),…,fp(sm))

已知点的协方差矩阵 Rij=R(θ,si,sj),i,j=1,2,…,m

未知点和已知点的协方差矩阵为 r(x)=[R(θ,s1,x),…,R(θ,sm,x)]T

我们希望用已知点的函数值来预测未知点的函数值，可以考虑如下的模型：

y^(x)=cTy

为了得到合适的系数，我们可以查看它的误差：

y^(x)−y(x)=cTy−y(x)=cT(Fβ+Z)−(f(x)Tβ+z)=cTZ−z+(FTc−f(x))Tβ

其中 Z=[z1,z2,…,zm]T 是已知点的误差,z是未知点的误差。 FTc 是已知点加权求和得到的未知点的基函数，应该和未知点的基函数相等，即

FTc(x)=f(x)

有了上述条件，现在我们就可以讨论预测的MSE了：

φ(x)=E(y^(x)−y(x))=E[(cTZ−z)2]=E(z2+cTZZTc−2cTZz)=σ2(1+cTRc−2cTr)

结合MSE和等式约束，可以构造拉格朗日乘子：

L(c,λ)=σ2(1+cTRc−2cTr)−λT(FTc−f)

∂L∂c=2σ2(Rc−r)−Fλ

所以可以得到如下方程：

[RFFT0][cλ^]=[rf]

λ^=−λ2σ2

解上述方程，得：

{λ^=(FTR−1F)−1(FTR−1r−f)c=R−1(r−Fλ^)

至此，推导中最重要的一步完成了，我们求出了c，就可以预测未知点：

y^(x)=cTY

4.3.3 通用形式

y^(x)=cTY

{λ^=(FTR−1F)−1(FTR−1r−f)c=R−1(r−Fλ^)

这个形式已经将预测模型求出，为了方便编程和理解，我们将进一步简化。将 λ^ 和c

带入表达式，得：

y^(x)=(r−Fλ^)TR−1Y=rTR−1Y−(FTR−1r−f)T(FTR−1F)−1FTR−1Y

假设 β∗=(FTR−1F)−1FTR−1Y ，上式可化简得：

y^(x)=rTR−1Y−(FTR−1r−f)Tβ∗=fTβ∗+rTR−1(Y−Fβ∗)=f(x)Tβ∗+r(x)Tγ∗

其中 Rγ∗=Y−Fβ∗

至此，我们已经推出了Kriging的表达式，任务完成！

4.3.4 重写MSE

将 和λ^和c 带入之前的MSE，得：

φ(x)=σ2(1+cT(Rc−2r))=σ2(1+(Fλ^−r)TR−1(Fλ^+r))=σ2(1+λ^TFTR−1Fλ^−rTR−1r)=σ2(1+uT(FTR−1F)−1u−rTR−1r)

其中， u=FTR−1r−f

得证！

4.3.5 求 σ

我们之前的推导求出了克里金模型的预测和方差表达式，表达式中的 σ 一直还不知道怎么求。接下来我们就详细推导 σ 应该等于多少。

我们的基本假设是目标函数可以用一系列基函数的加权求和得到，即

y=Fβ+e

所以数据点所对应的似然概率是

L=1(2πσ2)n/2|R|1/2exp[−(y−Fβ)TR−1(y−Fβ)2σ2]

对数似然概率为

ln(L)=−n2ln(2π)−n2ln(σ2)−12ln|R|−(y−Fβ)TR−1(y−Fβ)2σ2

为 了求得超参数 β,σ 的值，我们还是对上述似然函数求导，令梯度为0，由于 β=(β1,β2,…,βq)T 是一个向量。为了简化推导，我们假设 β 是一维向量，即只有一个基函数，我们对 β 求导，得

∂lnL∂β=−12σ2∂(y−Fβ)TR−1(y−Fβ)∂μ=0⟺∂(y−Fβ)TR−1(y−Fβ)∂β=0

(y−Fβ)TR−1(y−Fβ)=∑i=1m∑j=1m(yi−Fiβ)(yj−Fjβ)rij=∑i=1m∑j=1m(FiFjβ2−Fiyjβ−Fjyiβ+yiyj)rij

对上式求导，得：

2∑i=1m∑j=1mFirijFjβ−∑i=1m∑j=1mFirijyj−−∑j=1m∑i=1mFjrjiyi=0

2FTR−1Fβ=FTR−1Y+FTR−1Y

所以使得似然概率最大的 β 是

β∗=(FTR−1F)−1FTR−1Y

我们再对 σ 求导，将 β∗ 代入，得：

∂lnLσ=−m22σσ2−−22σ3(y−Fβ∗)TR−1(y−Fβ∗)=0

所以

σ2=1m(Y−Fβ∗)TR−1(Y−Fβ∗)

得证！

与形式一相同，在优化超参数时，使用的法则还是最大似然概率：

lnL≈−m2ln(σ2)−12ln(|R|)

4.3.6 DACE包的具体实现

尽管我们已经推导出了克里金模型的公式，在DACE包中的实际实现上时不是直接对上述矩阵求逆，而是加入了cholesky分解与QR分解来减少计算量。为了看懂MATLAB的代码，我们得到对最后的形式做一下变形。

1.求 β∗

首先对R做cholesky分解， C=chol(R),C=CT ，所以

β∗=(FT(CCT)F)−1FT(CCT)−1Y=((C−1F)T(C−1F))−1(C−1F)TC−1Y=(C−1F)−1C−1Y

取Ft=C\F, Yt=C\Y, [Q,G]=qr(Ft)，得

β∗=Ft−1Yt=(QG)−1Yt=G−1QTYt

2.求 γ∗

CCT=Y−Fβ∗CTγ∗=C−1Fβ∗=Yt−Ftβ∗=rho(define)

γ∗=C−Trho

3.求似然概率L

1m||rho||2=1m||Yt−Ftβ∗||2=1m||C−1Y−C−1Fβ∗||2=1m||C−1(Y−Fβ∗)||2=1m(Y−Fβ∗)TC−TC−1(Y−Fβ∗)=1m(Y−Fβ∗)TR−1(Y−Fβ∗)=σ2

L∼σ2∗|R|1m=σ2∗|C|2m

对照DACE中dacefit的源代码，可以很好地理解这一节的变换。

4.求预测值和方差的导数

形式二的好处在于将f(x)和r(x)分开了，因此预测值的导数很好求得：

∂y(x)∂x=Jf(x)β∗+Jr(x)Tγ∗

这里的 和Jf和Jr 是f和r的Jacobian矩阵

Jf(x)ij=∂fi∂xj(x),Jr(x)ij=∂R∂xj(θ,si,x)

 

s22(x)=σ2(1+uT(FTR−1F)−1u−rTR−1r)

同样，方差的导数同样是对x求导，得：

∂s2(x)∂x=2σ2(uT(FTR−1F)−1∂u∂x−rTR−1∂r∂x)=2σ2(uT(FTR−1F)−1(FTR−1∂f∂x−∂r∂x)−rTR−1∂r∂x)=2σ2[(uT(FTR−1F)−1FTR−1−rTR−1)∂r∂x−uT(FTR−1F)−1∂f∂x]

以上就是预测值与方差的导函数了。在DACE的具体实现中，仍然采用了QR分解等操作。

最后我们就介绍一下DACE中是如何减少 ∂s2(x)∂x 的计算量的。

rt=C−1r

u=FtTrt−f=FTC−TC−1r=FTR−1r

v=G−1u

Gv=G−Tu=G−TG−1u=(GGT)−1u

DACE只有在预测点数为1的时候才会给出方差的导数，此时G是一个数字，所以

GGT=GTG

又因为

FTR−1F=FT(CCT)−1F=FTC−TC−1F=(C−1F)T(C−1F)=(QG)T(QG)=GTQTQG=GTG=GGT

所以Gv可以重新写为

Gv=(FTR−1F)−1u

∂s2(x)∂x=2σ2[(FtGv−rt)TC−1dr−GvTdf]=2σ2[(C−1F(FTR−1F)−1u−C−1r)TC−1dr−((FTR−1F)−1u)Tdf]=2σ2[uT(FTR−1F)−1FTC−TC−1dr+rTC−TC−1dr−uT(FTR−1F)−1df]=2σ2[(uT(FTR−1F)−1FTR−1−rTR−1)∂r∂x−uT(FTR−1F)−1∂f∂x]

得证！