高斯过程回归【详细数学推导】

机器学习笔记

第一章 机器学习简介
第二章 感知机
第三章 支持向量机
第四章 朴素贝叶斯分类器
第五章 Logistic回归
第六章 线性回归和岭回归
第七章 多层感知机与反向传播【Python实例】
第八章 主成分分析【PCA降维】
第九章 隐马尔可夫模型
第十章 奇异值分解
第十一章 熵、交叉熵、KL散度
第十二章 什么是范数【向量范数、矩阵范数】
第十三章 熵、交叉熵、KL散度
第十四章 极大似然估计、最大后验估计、贝叶斯估计
第十五章 高斯过程回归模型

高斯过程回归

高斯过程回归（Gaussian Process Regression，简称GPR）是一种非参数的回归方法，它基于高斯过程（Gaussian Process）理论。在机器学习中，高斯过程是一种强大的工具，用于建模连续型的函数关系，特别适用于小样本的情况。在进行预测时，高斯过程回归会计算给定输入下目标函数的后验分布，这个后验分布也是一个高斯分布，它的均值给出了预测值，方差则表示了预测的不确定性。由于高斯过程回归是基于贝叶斯推断的，因此它能够提供预测结果的不确定性估计，这对于许多应用场景非常有用，尤其是在决策制定中需要考虑风险的情况下。

高斯过程回归的优点包括：

1. 能够灵活地处理非线性关系，适用于各种类型的函数逼近。
2. 能够提供预测的不确定性，有利于进行决策制定。
3. 不需要事先指定复杂的模型结构，减轻了参数调节的负担。

一、 一元高斯分布

若随机变量 $x$ 服从一个位置参数为 $\mu$ 、尺度参数为 $\sigma$ 的概率分布, 且其概率密度函数为 :
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$
则这个随机变量就称为正态随机变量, 正态随机变量服从的分布就称为正态分布(也称为高斯分布), 记作 $x\sim N\left(\mu ,{\sigma }^2\right)$。当 $\mu =0,\sigma =1$ 时, 正态分布就成为标准正态分布:
$$f(x)=\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{x^2}{2}\right)$$

一元高斯分布有很多很好的性质，比如两个高斯分布的线性组合仍然是高斯分布，而多个一元高斯分布的联合分布则可以导出多元高斯分布。

二、多元高斯分布

若一个多维随机向量具有同一元高斯分布类似的概率规律时, 称此随机向量遵从多元高斯(multivariate Gaussian)分布。多元高斯分布可以从一元高斯分布导出，对于 $D$ 维的随机变量量 $\mathbf{x}={\left(x_1,\dots ,x_D\right)}^T$ , 多元高斯分布形式为:
$$N(\mathbf{x}\mid \mu ,\mathbf{\Sigma })=\frac{1}{(2\pi {)}^{D/2}}\frac{1}{|\mathbf{\Sigma }{|}^{1/2}}\exp \left\{-\frac{1}{2}(\mathbf{x}-\mu {)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mu )\right\}$$
式中的 $\mu$ 是 $D$ 维的均值向量, $\Sigma$ 是 $D\times D$ 的协方差矩阵 (covariance matrix), $|\Sigma |$是$\Sigma$的行列式。多元正态分布也有很好的性质, 例如, 多元正态分布的边缘分布仍为正态分布，它经任何线性变换得到的随机向量仍为多维正态分布，它的线性组合为一元正态分布。在后面高斯过程回归的推导中会用到多元高斯分布的条件概率分布，所以这里我们给出多元高斯分布条件概率分布的定理：

定理（多元高斯分布条件概率分布）

设 $\mathbf{y}\sim \mathcal{N}(\mathbf{\mu },\Sigma )$, 将 $\mathbf{y},\mathbf{\mu }$ 和 $\Sigma$ 分割成如下形式：
$$\begin{array}{rl}& \mathbf{y}=\left[\begin{array}{l}{\mathbf{y}}_1\\ {\mathbf{y}}_2\end{array}\right]\begin{array}{l}{n}_1\\ n_2\end{array}\mathbf{\mu }=\left[\begin{array}{l}{\mathbf{\mu }}_1\\ {\mathbf{\mu }}_2\end{array}\right]\begin{array}{l}{n}_1\\ n_2\end{array}\\ & \Sigma =\left[\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right]\begin{array}{l}{n}_1\\ n_2\end{array}\end{array}$$
其中 $\mathbf{y},\mathbf{\mu }\in {\mathbb{R}}^n$ 且 $\Sigma \in {\mathbb{R}}^{n\times n}$。那么我们可以知道 ${\mathbf{y}}_1\sim \mathcal{N}\left({\mathbf{\mu }}_1,{\Sigma }_{11}\right)$ ，进一步可以推出条件概率分布：
$${\mathbf{y}}_2\mid {\mathbf{y}}_1\sim \mathcal{N}\left({\mathbf{\mu }}_2+{\Sigma }_{21}{\Sigma }_{11}^{-1}\left({\mathbf{y}}_1-{\mathbf{\mu }}_1\right),{\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{12}\right)$$

三、 高斯过程

高斯过程是概率论和数理统计中随机过程的一种，是一系列服从高斯分布的随机变量在一指数集内的组合。高斯过程中任意个随机变量的线性组合都服从多元高斯分布，每个有限维分布都是联合高斯分布。对一个任意集合 $\mathcal{X}$, 一个定义在$\mathcal{X}$上的高斯过程 $(\mathrm{GP})$ 是一个随机变量的集合 $(f(x),x\in \mathcal{X})$ ，使得对任意的 $n\in \mathbb{N}$ 且$x_1,\dots ,x_n\in \mathcal{X}$，满足$\left(f\left(x_1\right),\dots ,f\left(x_n\right)\right)$ 是一个多元高斯分布。由于高斯分布由均值向量和协方差矩阵指定，因此,高斯过程也完全由均值函数$\mu (x)$和协方差函数（或者叫核函数）$k(x,x^{\prime })$共同唯一确定，常表示为如下形式：
$$f(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(\mu (\mathbf{x}),k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)\right)$$
其中均值函数反映了函数在输入点x处的均值：
$$\mu (x)=\mathbb{E}(f(\mathbf{x}))$$
核函数$k(x,x^{\prime })$刻画了函数值在$x$和$x^{\prime }$处的关系：
$$k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=\mathbb{E}\left[(f(\mathbf{x})-\mu (\mathbf{x})\left(f\left({\mathbf{x}}^{\prime }\right)-\mu \left({\mathbf{x}}^{\prime }\right)\right)\right]$$
高斯过程的性质与其核函数有密切联系，不同的核函数可以赋予高斯过程不同的平滑性、各向同性、周期性和平稳性。核函数的选择要求满足Mercer定理(Mercer’s theorem)，即核函数在样本空间内的任意格拉姆矩阵为半正定矩阵。

四、高斯过程回归

我们知道在机器学习中，监督学习可以分为回归问题和分类问题。分类问题的输出是离散的类标签，而回归是关于连续量的预测。而本文，我们主要是讨论多任务高斯过程在回归问题中的应用，所以这里首先介绍一下高斯过程回归的基本原理。

高斯过程回归模型通常可以从权重空间视角和函数空间视角去推导，这里我们从函数空间的视角进行介绍。通过刚才对高斯过程的介绍，我们可以发现回归模型能用高斯过程来表示，比如贝叶斯线性回归模型$f(\mathbf{x})=\mathbf{\varphi }(\mathbf{x}{)}^{\mathrm{T}}\mathbf{w}$ ,其中$\mathbf{w}\sim \mathcal{N}\left(0,{\Sigma }_p\right)$、$\mathbf{\varphi }(\mathbf{x})$表示经过某种变换后的输入，于是可以得到：
$$\begin{array}{rl}\mu (\mathbf{x})& =\mathbb{E}[f(\mathbf{x})]=\mathbf{\varphi }(\mathbf{x}{)}^{\mathrm{T}}\mathbb{E}[\mathbf{w}]=0\\ k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)& =\mathbb{E}\left[f(\mathbf{x})f\left({\mathbf{x}}^{\prime }\right)\right]=\mathbf{\varphi }(\mathbf{x}{)}^{\mathrm{T}}\mathbb{E}\left[\mathbf{w}{\mathbf{w}}^{\mathrm{T}}\right]\mathbf{\varphi }\left({\mathbf{x}}^{\prime }\right)=\mathbf{\varphi }(\mathbf{x}{)}^{\mathrm{T}}{\Sigma }_p\mathbf{\varphi }\left({\mathbf{x}}^{\prime }\right)\end{array}$$
所以 $f(\mathbf{x})\sim \mathcal{G}\mathcal{P}\left(0,\mathbf{\varphi }(\mathbf{x}{)}^{\mathrm{T}}{\Sigma }_p\mathbf{\varphi }\left({\mathbf{x}}^{\prime }\right)\right)$。

进一步，我们可以考虑一个一般的带噪声的的回归模型：
$$y=f(\mathbf{x})+\epsilon f(\mathbf{x})\sim \mathcal{G}\mathcal{P}(\mu ,k)$$其中，$f(\mathbf{x})$是参数待定的高斯过程，噪声满足$\epsilon \sim \mathcal{N}\left(0,{\sigma }_n^2\right)$.现在我们假设数据集 D = { ( x i , y i ) ∣ i = 1 , … , n } \mathcal{D}=\left\{\left(\boldsymbol{x}_i, y_i\right) \mid i=1, \ldots, n\right\} D={(xi​,yi​)∣i=1,…,n}，其中 ${\mathbf{x}}_i\in {\mathbb{R}}^d,y_i\in \mathbb{R}$。可以表示为矩阵形式$\mathcal{D}=(\mathrm{X},\mathbf{y})$，其中$\mathrm{X}\in {\mathbb{R}}^{n\times d},\mathbf{y}\in {\mathbb{R}}^n$。当观测点满足回归模型时，按照高斯过程的定义,这些点的联合分布 $\left[f\left({\mathbf{x}}_1\right),\dots ,f\left({\mathbf{x}}_n\right)\right]$ 需要满足一个多维高斯分布, 即:
$${\left[f\left({\mathbf{x}}_1\right),f\left({\mathbf{x}}_2\right),\dots ,f\left({\mathbf{x}}_n\right)\right]}^{\mathrm{T}}\sim \mathcal{N}(\mathbf{\mu },K)$$这里 $\mathbf{\mu }={\left[\mu \left({\mathbf{x}}_1\right),\dots ,\mu \left({\mathbf{x}}_n\right)\right]}^{\mathrm{T}}$ 是均值向量, $K$ 是 $n\times n$ 的矩阵, 其中第 $(i,j)$ 个元素是 $K_{ij}=k\left({\mathbf{x}}_i,{\mathbf{x}}_j\right)$.

为了预测在新的点处的函数值 $f_{\ast }=f(Z)$, 其中 $Z={\left[{\mathbf{z}}_1,\cdots ,{\mathbf{z}}_m\right]}^{\mathrm{T}}$, 根据高斯分布的性质, 训练点和预测点的联合分布仍然为高斯分布，即:
$$\left[\begin{array}{c}\mathbf{y}\\ f_{\ast }\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}\mathbf{\mu }(X)\\ \mathbf{\mu }(Z)\end{array}\right],\left[\begin{array}{cc}K(X,X)+{\sigma }_n^2\mathbf{I}& K(Z,X{)}^{\mathrm{T}}\\ K(Z,X)& K(Z,Z)\end{array}\right]\right)$$
其中:

• $\mathbf{\mu }(X)=\mathbf{\mu },\mathbf{\mu }(Z)={\left[\mu \left({\mathbf{z}}_1\right),\dots ,\mu \left({\mathbf{z}}_m\right)\right]}^{\mathrm{T}},K(X,X)=K$;
• $K(Z,X)$是个 $m\times n$ 的矩阵, 其中第 $(i,j)$ 个元素 $[K(Z,X){]}_{ij}=k\left({\mathbf{z}}_i,{\mathbf{x}}_j\right)$ ；
• $K(Z,Z)$是个 $m\times m$ 的 矩阵, 其中第 $(i,j)$ 个元素 $[K(Z,Z){]}_{ij}=k\left({\mathbf{z}}_i,{\mathbf{z}}_j\right)$ .
  最后利用高斯分布的条件分布性质, 我们可以得到关于预测值的条件概率分布:
  $$p\left(f_{\ast }\mid X,\mathbf{y},Z\right)=\mathcal{N}(\hat{\mathbf{\mu }},\hat{\Sigma })$$其中：
  $$\begin{array}{rl}& \hat{\mathbf{\mu }}=K(Z,X){\left(K(X,X)+{\sigma }_n^2\mathbf{I}\right)}^{-1}(\mathbf{y}-\mathbf{\mu }(X))+\mathbf{\mu }(Z)\\ & \hat{\Sigma }=K(Z,Z)-K(Z,X){\left(K(X,X)+{\sigma }_n^2\mathbf{I}\right)}^{-1}K(Z,X{)}^T\end{array}$$如果将预测值的噪声考虑进来，条件概率分布如下：
  $$p\left({\mathbf{y}}_{\ast }\mid X,\mathbf{y},Z\right)=\mathcal{N}\left(\hat{\mathbf{\mu }},\hat{\Sigma }+{\sigma }_n^2\mathbf{I}\right)$$在实际应用中，我们通常将令均值函数$\mu (x)=0$(本文如无特殊说明，也都采用值为0的均值函数)，那么预测的均值函数和预测的协方差函数的结果可以更为简单的表示为:
  $$\begin{array}{rl}\hat{\mathbf{\mu }}& =K(Z,X)(K(X,X)+{\sigma }_n^2\mathbf{I}{)}^{-1}\mathbf{y}\\ \hat{\Sigma }& =K(Z,Z)-K(Z,X)(K(X,X)+{\sigma }_n^2\mathbf{I}{)}^{-1}K(Z,X{)}^{\mathrm{T}}\end{array}$$
  至此，我们便完成了从函数空间视角出发的高斯过程回归模型的推导。观察上式，我们可以发现一些有趣的性质：

1. 首先来看均值$\hat{\mathbf{\mu }}$，由于测试数据集一共有m个点，因而 $\hat{\mathbf{\mu }}$理应是 $m\times 1$，而对应等式右边的 $\mathbf{y}$应当是 n×1 ，而$\hat{\mathbf{\mu }}$等式右边除了$\mathbf{y}$的其他部分理应为 m×n的，所以预测均值是观测点$\mathbf{y}$的线性组合。
2. 再来看协方差$\hat{\mathbf{\Sigma }}$，等式右边的第一部分是我们的先验的协方差，减掉的后面的那一项实际上表示了观测到数据后函数分布不确定性的减少，如果第二项非常接近于 0，说明观测数据后我们的不确定性几乎不变，反之如果第二项非常大，则说明不确定性降低了很多。

如果我们再换一种视角来看均值$\hat{\mathbf{\mu }}$，将${\left(K(X,X)+{\sigma }_n^2\mathbf{I}\right)}^{-1}\mathbf{y}$ 看做整体 $\alpha$, 则
$${\hat{\mathbf{\mu }}}_i=\sum _{j=1}^n{\alpha }_{j}k\left({\mathbf{x}}_j,{\mathbf{z}}_i\right)$$最终结果表示，预测值${\hat{\mathbf{\mu }}}_i$可以看做是预测点和观测点之间核函数的线性组合，这个方程表明高斯过程回归等价于使用核函数$k(x,x^{\prime })$将输入投影到特征空间的线性回归模型。为了得到预测点的值，每个输出$\widehat{{\mathbf{\mu }}_i}$被加权为其相关的输入$x_j$与被预测点$z_i$的相似度，这是由核函数诱导的相似度度量。这个加权和表明，高斯过程回归模型取决于选定的核函数和到目前为止观察到的数据，而不需要指定具体的模型函数形式，这就是为什么高斯过程回归是一种非参数方法。

用高斯过程对回归问题进行建模有很多优势。首先，GP利用有限的训练数据点计算核函数，建模复杂度与 $x$的维数无关 ，使得我们能够处理更高维的输入数据，理论上无限维都行。其次，刚刚提到，高斯过程是一种非参数的方法，对回归问题进行建模时，不需要指定模型的形式，这意味着高斯过程可以对任意形式的函数进行建模拟合。

五、核函数

我们看到$k(x,x^{\prime })$对高斯过程回归模型至关重要，可以说核函数是一个高斯过程的核心，核函数决定了一个高斯过程的性质。

定义（核函数）

在$\mathcal{X}$上的半正定核函数$k$，需要满足如下条件 $k:\mathcal{X}\times \mathcal{X}\mapsto \mathbb{R}$, $\forall n\in \mathbb{N},\forall x_1,\dots ,x_n\in \mathcal{X}$, s.t.矩阵 $C$ 都是半正定的, 其中 $C_{ij}=k\left(x_i,x_j\right)$

比如，当 $\mathcal{X}={\mathbb{R}}^d,k(\mathbf{x},{\mathbf{x}}^{\prime })={\mathbf{x}}^{\mathrm{T}}{\mathbf{x}}^{\prime }$时，可以得到 $C=\mathbf{x}{\mathbf{x}}^{\mathbf{T}}$. 设 $\mathbf{a}\in {\mathbb{R}}^n$, 则 ${\mathbf{a}}^{\mathrm{T}}\mathbf{C}\mathbf{a}={\mathbf{a}}^{\mathrm{T}}\mathbf{x}{\mathbf{x}}^{\mathrm{T}}\mathbf{a}={\left(\mathbf{a}{\mathbf{x}}^{\mathrm{T}}\right)}^2\ge 0$. 所以，二元函数$k(\mathbf{x},\mathbf{y})={\mathbf{x}}^{\mathrm{T}}\mathbf{y}$ 是一个半正定的核函数。

在高斯过程中，核函数用于生成一个协方差矩阵来衡量任意两个点之间的“距离”，一个合理的假设通常是两点之间的相关性随两点之间的距离而衰减。 不同的核函数有不同的衡量方法，得到的高斯过程的性质也不一样。最常用的一个核函数为高斯核函数，也称为平方指数函数(Squared Exponential,SE)，其基本形式如下：
$$k(x_i，x_j)=s_f^{2}exp(-\frac{\Vert x_i-x_j{\Vert }_2^2}{2l^2})$$其中$s_f$和$l$是高斯核的超参数，$s_f$是信号方差(signal variance)，控制输出的振幅。$l$是长度尺度(length scale)，控制输入方向的振荡频率，下图展示了两个超参数对高斯过程函数性质的影响。

左边图展示了，当$s_f$固定时，长度尺度$l$对高斯过程的影响，可以看到不同的$l$导致高斯过程曲线在水平轴上具有不同的振荡频率。右边图展示了，当$l$固定时,信号方差$s_f$对高斯过程的影响，可以看到不同的$s_f$导致高斯过程曲线在纵轴上有不同的振幅。

高斯核函数提供了一个富有表现力的核来建模平滑、平稳的函数，超参数$l$和$s_f$可以控制点与点之间的先验相关性，从而控制高斯过程的预测性能。

六、超参数的优化

通过前面的介绍，我们已经知道了如何使用给定的核和零均值函数构造高斯过程回归模型。只要从数据中学习所有待定超参数 θ = { θ 1 , θ 2 , ⋯   } \boldsymbol{\theta}=\{\theta_1,\theta_2,\cdots\} θ={θ1​,θ2​,⋯}，就可以得到预测均值和方差。

由于高斯过程回归模型中超参数的后验分布不易获得，因此在实际应用中不经常使用超参数的全贝叶斯推理，通常的做法是通过最大化边际似然来获得超参数的点估计。所以可以采用最大似然估计(Maximum Likelihood Estimation，MLE)，通过梯度下降等优化方法进行优化。

高斯过程模型的似然函数可以表示为:
$$p(\mathbf{y}\mid X,\theta )=\int p(\mathbf{y}\mid \mathbf{f},X,\theta )p(\mathbf{f}\mid X,\theta )\mathrm{d}\mathbf{f}$$
在高斯过程回归模型中，先验和似然函数都是高斯分布：
$$\begin{array}{rl}p(f\mid X,\mathbf{\theta })& =\mathcal{N}(0,K),\\ p(\mathbf{y}\mid f,X,\mathbf{\theta })& =\mathcal{N}\left(f,{\sigma }_n^2\mathrm{I}\right).\end{array}$$于是我们知道边缘似然函数也是服从高斯分布：
$$p(\mathbf{y}\mid X,\mathbf{\theta })=\int \mathcal{N}\left(f,{\sigma }_n^2\mathrm{I}\right)\mathcal{N}(0,K)\mathrm{d}f=\mathcal{N}\left(0,K+{\sigma }_n^2\mathrm{I}\right)=\mathcal{N}\left(0,{\Sigma }_{\theta }\right),$$
其中 ${\Sigma }_{\theta }=K_{\theta }+{\sigma }_n^2\mathbf{I}=K+{\sigma }_n^2\mathbf{I}$ ， $\mathbf{\theta }$ 包含于协方差函数 $K$中。而更常用的是负对数边缘似然函数(Negative Log Marginal Likelihood,NLML)：
$$\mathcal{L}(\mathbf{\theta })=-\log p(\mathbf{y}\mid X,\mathbf{\theta })=\frac{1}{2}{\mathbf{y}}^{\mathrm{T}}{\Sigma }_{\theta }^{-1}\mathbf{y}+\frac{1}{2}\log |{\Sigma }_{\theta }|+\frac{n}{2}\log 2\pi$$
得到负对数似然函数$\mathcal{L}$后，剩下的就是优化问题了，可以对上式关于$\theta$求偏导：
$$\frac{\partial \mathcal{L}}{\partial {\theta }_i}=\frac{1}{2}\mathrm{tr}\left({\Sigma }_{\theta }^{-1}\frac{\partial {\Sigma }_{\theta }}{\partial {\theta }_i}\right)-\frac{1}{2}{\mathbf{y}}^{\mathrm{T}}{\Sigma }_{\theta }^{-1}\frac{\partial {\Sigma }_{\theta }}{\partial {\theta }_i}{\Sigma }_{\theta }^{-1}\mathbf{y}.$$

于是，就可以利用梯度下降之类的方法求得超参数$\mathbf{\theta }$。不过值得注意的是，MLE方法有一些局限性。对于许多核函数，边际似然函数相对于超参数不是凸的，因此优化算法可能收敛到局部最优点，而不是全局最优点。因此，通过最大似然估计优化超参数，得到的GPR性能取决于优化算法的初始值 ，所以如何对模型参数的初始化十分关键。

参考资料

1. 高斯过程的简洁推导