高斯过程回归（输出学习法！）

RL实验跑多了，理论就基本忘光光了，只剩下intuitive的视角，本想着用贝叶斯优化来调超参亦或者指导RL中的探索，但第一个门槛是对高斯过程的直觉理解一点都不深入。看着网上许多高斯过程回归的竟然看不太懂，然后还是自己写写吧。主要浓缩自2006年的参考书目Gaussian Process for Machine Learning的Chapter 2以及一篇英文blog visual gaussian process

概述

首先一个监督问题有三个对象： 数据$X$、模型(结构与参数$W$)、标签$Y$，一般有数据集 { X , Y } \{X,Y\} {X,Y}，然后我们需要选择模型结构（线性linear、平滑smoothing、凸凹等特性）以及学习算法(MLE、MAP)来得到参数$W$.

因此解决监督问题一般有两个方向：

1. 限制模型结构的选择（比如选线性、平滑性好的模型结构)
2. 对所有可能的模型结构设定一个prior（希望能从data中学习适应自身结构的模型结构，如适合的给high probability）

对于第一个方向，诸如L1或L2等regularization的方式都属于此类，但因为constraint的存在，容易出现overfitting的现象，需要richness of function class来寻找更适合Data自身结构的模型
对于第二个方向，所有可能的模型结构(function class)是无穷多个的，阿这？咋搞？这时候就搬来了Gaussian Process。这也就是说高斯过程是由定义在连续域上的服从高斯分布的无穷维随机变量组成的。

多元高斯分布的定义

---

一个n维随机向量$\mathbf{X}$服从多元高斯分布，均值向量$\mu$决定了分布的中心，协方差矩阵$\Sigma$决定了分布的形状且反映了各维度变量之间的关系：

$$\begin{gathered} \mathbf{X}=\left[\begin{array}{c}{X}_1\\ X_2\\ ⋮\\ X_n\end{array}\right]\sim \mathcal{N}(\mu ,\Sigma ) \\ where{X}_i\sim \mathcal{N}({\mu }_i,{\Sigma }_i) \\ \Sigma =Cov(X_i,X_j)=E[(X_i-{\mu }_i)(X_j-{\mu }_j)] \end{gathered}$$

多元高斯分布具有下述的性质：

边际分布：$$X\sim \mathcal{N}({\mu }_X,{\Sigma }_X),Y\sim \mathcal{N}({\mu }_Y,{\Sigma }_Y)$$

联合分布为：

$$P_{X,Y}=P(x,y)=\left[\begin{array}{c}X\\ Y\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mu }_X\\ {\mu }_Y\end{array}\right],\left[\begin{array}{cc}{\Sigma }_{XX}& {\Sigma }_{XY}\\ {\Sigma }_{YX}& {\Sigma }_{YY}\end{array}\right]\right)$$

条件分布：

$$\begin{array}{l}X\mid Y\sim \mathcal{N}\left({\mu }_X+{\Sigma }_{XY}{\Sigma }_{YY}^{-1}\left(Y-{\mu }_Y\right),{\Sigma }_{XX}-{\Sigma }_{XY}{\Sigma }_{YY}^{-1}{\Sigma }_{YX}\right)\\ Y\mid X\sim \mathcal{N}\left({\mu }_Y+{\Sigma }_{YX}{\Sigma }_{XX}^{-1}\left(X-{\mu }_X\right),{\Sigma }_{YY}-{\Sigma }_{YX}{\Sigma }_{XX}^{-1}{\Sigma }_{XY}\right)\end{array}$$

需要明确的是，一般来说这里的随机变量是定义在维度dimension上的，即一般是属性attribute（如影响商品价格的attribute有质量$X_1$，供求比$X_2$，成本$X_3$等），所以价格服从有限多元的联合高斯分布。

高斯过程

---

定义如下：
一个高斯过程是无限多个服从高斯分布随机变量的集合，从该集合中随意选取的有限个随机变量服从自身的联合高斯分布。
假设取三个随机变量$i,j,k$即$X_1,...X_n,n\to \infty ,\forall i,j,k\in n,P_{X_i,X_j,X_k}\sim \mathcal{N}({\mu }_{ijk},{\Sigma }_{ijk})$，因此从数学意义上，无限个随机变量 { X n } n = 1 n → ∞ \{X_n\}_{n=1}^{n\rightarrow\infty} {Xn​}n=1n→∞​可以看作是多元高斯分布的维度dim拉长到无限，从函数的角度用$f(x)$进行代替，描述为：

$$\begin{array}{rl}f(\mathbf{x})& \sim \mathcal{G}\mathcal{P}\left(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime })\right)\\ wherem(\mathbf{x})& =\mathbb{E}[f(\mathbf{x})]\\ k(\mathbf{x},{\mathbf{x}}^{\prime })& =\mathbb{E}[(f(\mathbf{x})-m(\mathbf{x}))(f({\mathbf{x}}^{\prime })-m({\mathbf{x}}^{\prime }))]\end{array}\text{\hspace{1em}(*)}$$

尽管数学意义上是有限维度到无限维度上的拓展，但实际使用的时候，随机变量的含义就发生了改变。使用多元高斯分布时，一个随机变量是一个维度上的高斯分布，因此可以将一维度对应一个attribute，维度是离散的；但高斯过程是定义在连续域上的，一个随机变量是连续域上的一个值，维度是连续的。最常见的情况，随机变量的index所在的连续域是时间。

这里初看定义，还是有点懵，直觉有阻碍，下面从两个角度阐释高斯过程。

1. Weight-Space View即从模型的参数$W$的角度出发，看待高斯过程
2. Function-space View即把高斯过程看作是函数的分布，从函数空间的角度进行推断

一、Weight-space View

1.1 标准贝叶斯线性模型

此处回顾一下标准的多元线性回归模型 + 贝叶斯

• 有$n$个样本的数据集 D = { ( x i , y i ) } i = 1 n \mathcal D=\{(\mathbf{x_i},y_i)\}_{i=1}^n D={(xi​,yi​)}i=1n​，其中$x_i\in {\mathbb{R}}^D,y_i\in \mathbb{R}$，输入矩阵形式为 $X\in {\mathbb{R}}^{D\times n}$，输出targets向量为$\mathbf{y}\in {\mathbb{R}}^n$.
  这里输入矩阵的表示，使得一个样本是一列.

• 选择标准的线性模型（模型结构为线性，并受取决于样本的高斯噪声扰动）
  $$f(\mathbf{x})={\mathbf{x}}^T\mathbf{w},y=f(\mathbf{x})+ϵ,ϵ\sim \mathcal{N}(0,{\sigma }_n^2)$$

  因此其参数对象为$\mathbf{w}\in {\mathbb{R}}^D$，于是这个监督问题便被建模起来了，从$X,Y$中学习线性结构模型中的参数$W$即$p(\mathbf{w}|X,\mathbf{y})$

• 贝叶斯线性模型 ( 请自行对应posterior, likelihood, prior )
  $$\text{posterior}=\frac{\text{likelihood}\times \text{prior}}{\text{marginal likelihood}}⟺p(\mathbf{w}\mid \mathbf{y},X)=\frac{p(\mathbf{y}|X,\mathbf{w})p(\mathbf{w})}{p(\mathbf{y}\mid X)}$$

  频率派角度认为参数w是一个数，而贝叶斯派认为参数w服从一个分布，此处为贝叶斯学派的角度。于是参数分布通过线性结构(function class)决定了模型。MAP与MLE的差别仅在于prior，此处指定参数分布的prior服从均值为0向量，协方差矩阵为${\Sigma }_p$的多元高斯分布：
  $$\mathbf{w}\sim \mathcal{N}(0,{\Sigma }_p)$$

  likelihood为（假设样本i.i.d)：
  $$\begin{array}{rl}p(\mathbf{y}\mid X,\mathbf{w})& =\prod _{i=1}^{n}p\left(y_i\mid {\mathbf{x}}_i,\mathbf{w}\right)=\prod _{i=1}^n\frac{1}{\sqrt{2\pi }{\sigma }_n}\exp \left(-\frac{{\left(y_i-{\mathbf{x}}_i^{\top }\mathbf{w}\right)}^2}{2{\sigma }_n^2}\right)\\ & =\frac{1}{{\left(2\pi {\sigma }_n^2\right)}^{n/2}}\exp \left(-\frac{1}{2{\sigma }_n^2}{\left|\mathbf{y}-X^{\top }\mathbf{w}\right|}^2\right)=\mathcal{N}\left(X^{\top }\mathbf{w},{\sigma }_n^{2}I\right)\end{array}$$

  posterior为：

  $$\begin{gathered} \begin{array}{rl}p(\mathbf{w}\mid X,\mathbf{y})& \propto \exp \left(-\frac{1}{2{\sigma }_n^2}{\left(\mathbf{y}-X^{\top }\mathbf{w}\right)}^{\top }\left(\mathbf{y}-X^{\top }\mathbf{w}\right)\right)\exp \left(-\frac{1}{2}{\mathbf{w}}^{\top }{\Sigma }_p^{-1}\mathbf{w}\right)\\ & \propto \exp \left(-\frac{1}{2}(\mathbf{w}-\overline{\mathbf{w}}{)}^{\top }\left(\frac{1}{{\sigma }_n^2}X{X}^{\top }+{\Sigma }_p^{-1}\right)(\mathbf{w}-\overline{\mathbf{w}})\right)\end{array} \\ where\overline{\mathbf{w}}={\sigma }_n^{-2}({\sigma }_n^{-2}X{X}^T+{\Sigma }_p^{-1}{)}^{-1}X\mathbf{y} \end{gathered}$$

  令 $\mathbf{A}={\sigma }_n^{-2}X{X}^T+{\Sigma }_p^{-1}$有：
  $$p(\mathbf{w}|X,\mathbf{y})\sim \mathcal{N}(\overline{\mathbf{w}}=\frac{1}{{\sigma }_n^2}{A}^{-1}X\mathbf{y},A^{-1})$$

  因此通过数据集中的$X,\mathbf{y}$，指定模型$f(x)$的输出噪声为$ϵ\sim \mathcal{N}(0,{\sigma }_n^2)$，初始化模型参数的prior为$\mathbf{w}\sim \mathcal{N}(0,{\Sigma }_p)$，便可得到参数的分布$p(\mathbf{w}|X,\mathbf{y})$.

• 贝叶斯推断：给定测试样本$x_{\ast }$，推断其输出$f(x_{\ast })$，简记为$f_{\star }$

  $$\begin{array}{rl}p\left(f_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\right)& =\int p\left(f_{\ast }\mid {\mathbf{x}}_{\ast },\mathbf{w}\right)p(\mathbf{w}\mid X,\mathbf{y})d\mathbf{w}\\ & =\mathcal{N}\left(\frac{1}{{\sigma }_n^2}{\mathbf{x}}_{\ast }^{\top }{A}^{-1}X\mathbf{y},{\mathbf{x}}_{\ast }^{\top }{A}^{-1}{\mathbf{x}}_{\ast }\right)\end{array}$$

根据下图总结一下：

1. 具体化问题为一维的$x$，两个参数$w_1,w_2$的线性模型$f(x)=w_1+w_{2}x$，noise为$ϵ\sim \mathcal{N}(0,I)$
2. 确定prior为$w\sim \mathcal{N}(0,\mathbf{I})$，对应下图(a)
3. 确定数据集$X,\mathbf{y}$为图(b)中的三个带星号的样本点
4. 于是likelihood和posterior都有了$p(\mathbf{y}\mid X,\mathbf{w})=\mathcal{N}\left(X^{\top }\mathbf{w},{\sigma }_n^{2}I\right)$、$p(\mathbf{w}|X,\mathbf{y})\sim \mathcal{N}(\overline{\mathbf{w}}=\frac{1}{{\sigma }_n^2}{A}^{-1}X\mathbf{y},A^{-1})$，分别对应（c)(d)
5. 然后经过贝叶斯推断公式$p\left(f_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\right)=\mathcal{N}\left(\frac{1}{{\sigma }_n^2}{\mathbf{x}}_{\ast }^{\top }{A}^{-1}X\mathbf{y},{\mathbf{x}}_{\ast }^{\top }{A}^{-1}{\mathbf{x}}_{\ast }\right)$，在(b)中画出mean plot对应实线，加出一两个标准差对应虚线。

按逻辑抽象总结一波：

1. 将一定数量的样本组成$X$，对应的label组成$\mathbf{y}$
2. 将它们丢进设定好prior与model noise的标准线性模型中
3. 根据MAP准则得到模型参数$W$服从的分布，即posterior
4. 给一个测试样本$x_{\star }$，模型通过贝叶斯推断给label服从的分布$f_{\star }$
5. 从label分布中取mean值，or sample给出estimate

1.2 Feature Space的贝叶斯线性模型

上述之所以为标准线性模型，是因为定义模型时采用的是直接的输出：
$$f(\mathbf{x})={\mathbf{x}}^{T}w$$

该线性指的是与参数的线性，将$D$维输入向量$\mathbf{x}$投射到高维空间$\varphi (x)$，认为是特征空间Feature Space，在高维空间上建模与参数的线性模型，对原本输入其中第一维度$x_1$投射形式可以选多项式方式，扩展到$L$维。
$$\varphi (x_1)=(1,x,x^2,...,x^{L-1}{)}^T$$

无论哪种方式，$\varphi$的作用就是把$D$维输入向量$\mathbf{x}\in {\mathbb{R}}^D$中的每一维度亦不同的方式扩展成高维$N$维$\varphi (x)\in {\mathbb{R}}^N$.

贝叶斯线性模型为：
$$f(\mathbf{x})=\varphi (\mathbf{x}{)}^{T}w$$

只需要把1.1节贝叶斯标准线性模型中的推断形式换一下即可。
原来为：
$$\begin{array}{rl}p\left(f_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\right)& =\int p\left(f_{\ast }\mid {\mathbf{x}}_{\ast },\mathbf{w}\right)p(\mathbf{w}\mid X,\mathbf{y})d\mathbf{w}\\ & =\mathcal{N}\left(\frac{1}{{\sigma }_n^2}{\mathbf{x}}_{\ast }^{\top }{A}^{-1}X\mathbf{y},{\mathbf{x}}_{\ast }^{\top }{A}^{-1}{\mathbf{x}}_{\ast }\right)\end{array}\text{\hspace{1em}(1)}$$

将$X$换成$\mathbf{\Phi }(X)$和${\mathbf{x}}_{\ast }$换成$\mathbf{\varphi }\left({\mathbf{x}}_{\ast }\right)$，于是(1)变成:
$$f_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\sim \mathcal{N}\left(\frac{1}{{\sigma }_n^2}\mathbf{\varphi }{\left({\mathbf{x}}_{\ast }\right)}^{\top }{A}^{-1}\mathbf{\Phi }(X)\mathbf{y},\mathbf{\varphi }{\left({\mathbf{x}}_{\ast }\right)}^{\top }{A}^{-1}\mathbf{\varphi }\left({\mathbf{x}}_{\ast }\right)\right)\text{\hspace{1em}(2)}$$

其中$A={\sigma }_n^{-2}\mathbf{\Phi }(X)\mathbf{\Phi }(X{)}^T+{\Sigma }_p^{-1}$，${\Sigma }_p^{-1}$是来自prior的。

老生常谈的一点是，(2)需要对高维的矩阵$A$是$N\times N$求逆，需要很大的computation，需要引入kernel trick，怎么引入呢？

将$\mathbf{\Phi }(X)$,$\mathbf{\varphi }\left({\mathbf{x}}_{\ast }\right)$即简记为$\Phi ,{\varphi }_{\ast }$，于是$A={\sigma }_n^{-2}\Phi {\Phi }^{\top }+{\Sigma }_p^{-1}$，定义$K={\Phi }^{\top }{\Sigma }_p\Phi$，(2)式恒等变换成：
$$\begin{array}{rl}{f}_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\sim \mathcal{N}(& {\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi {\left(K+{\sigma }_n^{2}I\right)}^{-1}\mathbf{y},\\ & {\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast }-{\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi {\left(K+{\sigma }_n^{2}I\right)}^{-1}{\Phi }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast })\end{array}\text{\hspace{1em}(3)}$$

此处变换正说明了该算法适用范围是样本数量比较少，但样本维度特别高的情况！

该变换成功将需要求逆的矩阵从$N\times N$变成$n\times n$（n为样本数量），观察(3)式，有${\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi ,{\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast },{\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi ,{\Phi }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast }$这些仍需要计算高维inner product的项，将$n\times n$协方差矩阵${\Sigma }_p$进行SVD分解有${\Sigma }_p=UD{U}^{\top }$，${\Sigma }_p^{1/2}=U{D}^{1/2}{U}^{\top }$,因此定义kernel为：

$$\begin{array}{rl}k(\mathbf{x},{\mathbf{x}}_{\ast })& ={\Phi }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast }\\ & =({\Sigma }_p^{1/2}\Phi {)}^{\top }({\Sigma }_p^{1/2}{\mathbf{\varphi }}_{\ast })\\ & =\psi (\mathbf{x}{)}^{\top }\psi ({\mathbf{x}}_{\ast })\end{array}$$

所以(3)中分布的参数就可以全用这个kernel进行计算了。kernel trick确实是为了简化高维计算才引入的，其也塑造了原始输入映射到高维空间的特征构造方式。所以选择了一个kernel，就意味着特征映射到高维空间的方式确定了。

二、Function-space View

2.1 高斯过程回顾GPR前情提要

高斯过程的定义：高斯过程是随机变量的集合，随机变量的任何有限数量都服从联合高斯分布

$$\begin{array}{rl}f(\mathbf{x})& \sim \mathcal{G}\mathcal{P}\left(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime })\right)\\ wherem(\mathbf{x})& =\mathbb{E}[f(\mathbf{x})]\\ k(\mathbf{x},{\mathbf{x}}^{\prime })& =\mathbb{E}[(f(\mathbf{x})-m(\mathbf{x}))(f({\mathbf{x}}^{\prime })-m({\mathbf{x}}^{\prime }))]\end{array}\text{\hspace{1em}(*)}$$

Weight-Space View主要是想说明贝叶斯线性回归是高斯过程的特例。

将$f(\mathbf{x})=\mathbf{\varphi }(\mathbf{x}{)}^{\top }\mathbf{w},\mathbf{w}\sim \mathcal{N}(0,{\Sigma }_p)$代入高斯过程的定义(*)中：
$$\begin{array}{rl}\mathbb{E}[f(\mathbf{x})]& =\mathbf{\varphi }(\mathbf{x}{)}^{\top }\mathbb{E}[\mathbf{w}]=0\\ \mathbb{E}\left[f(\mathbf{x})f\left({\mathbf{x}}^{\prime }\right)\right]& =\mathbf{\varphi }(\mathbf{x}{)}^{\top }\mathbb{E}\left[\mathbf{w}{\mathbf{w}}^{\top }\right]\mathbf{\varphi }\left({\mathbf{x}}^{\prime }\right)=\mathbf{\varphi }(\mathbf{x}{)}^{\top }{\Sigma }_p\mathbf{\varphi }\left({\mathbf{x}}^{\prime }\right).\end{array}$$

所以贝叶斯线性回归是mean function即$m(\mathbf{x})=0$，covariance function为kernel的特例。
这从直觉上怎么理解？

1. mean function为0，意味着每一个样本$\mathbf{x}$通过模型$f$的输出$f(\mathbf{x})$期望为0，在没有特别先验的情况下，均可设置所有可能的函数在某个样本的值均值为0（这里的期望的对象是所有函数，在一个样本点$x_0$下，$f_1(x_0),f_2(x_0),...,f_n(x_0)$加起来为0）
2. covariance function取决于kernel，意味着模型输出$f(\mathbf{x})$与回归值$\mathbf{y}$的不同之处在于该样本$\mathbf{x}$与其它样本${\mathbf{x}}^{\prime }$的距离$\mathbf{\varphi }(\mathbf{x}{)}^{\top }{\Sigma }_p\mathbf{\varphi }\left({\mathbf{x}}^{\prime }\right)$，距离越远，模型输出$f(\mathbf{x})$与回归值$\mathbf{y}$的扰动越大。将样本的距离与输出空间的函数分布扰动建立了联系（符合离其它样本${\mathbf{x}}^{\prime }$越远的样本$\mathbf{x}$，其输出的$f(\mathbf{x})$应该更远离的直觉）

所以从函数的角度上看：

• mean function就是函数$f(x)$分布的prior
• covariance functoon就是prior中关于函数分布的平滑特性smoothness

2.2 指定covariance function的GPR

既然mean function是函数分布的prior，在一无所知的情况下一般指定为0，于是关键是选择一个特定的covariance function，此处选择指定如下：

$$\begin{array}{rl}\mathrm{cov}\left(f\left(\mathbf{x}\right),f\left({\mathbf{x}}^{\prime }\right)\right)& =\mathbb{E}[(f(\mathbf{x})-m(\mathbf{x}))(f({\mathbf{x}}^{\prime })-m({\mathbf{x}}^{\prime }))]\\ & =k\left(\mathbf{x},\mathbf{x}\right)=\exp \left(-\frac{1}{2}{\left|\mathbf{x}-{\mathbf{x}}^{\prime }\right|}^2\right)\\ & wherem(\mathbf{x})=m({\mathbf{x}}^{\prime })=0\end{array}$$

这covariance function是两个输入样本之间距离的函数，两个样本越近说明它们对彼此的影响越大即方差越大，越远影响越小从而方差越小，还可以引入一个length scale即$\ell$来衡量距离的程度对输出的影响。
$$\left|\mathbf{x}-{\mathbf{x}}^{\prime }\right|/\ell$$

所以总结一波，问题如下：$(x\in \mathbb{R})$
$$\begin{array}{rl}f(\mathbf{x})& \sim \mathcal{G}\mathcal{P}\left(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime })\right)\\ & wherem(\mathbf{x})=0\\ & k(\mathbf{x},{\mathbf{x}}^{\prime })=\exp \left(-\frac{1}{2}{\left|x-x^{\prime }\right|}^2\right)\end{array}$$

下面的测试点只是为了画出函数！请重点理解！

1. 在prior的情况下，从$[-5,5]$中选择一堆测试点$x_1,...,x_{n_{\ast }}$拼成矩阵$X_{\ast }$，然后逐个计算$k(x_p,x_q),\forall p,q\in n_{\ast }$拼成$n_{\ast }\times n_{\ast }$的矩阵$K(X_{\ast },X_{\ast })$，形成$$\mathbf{f}(X_{\ast })\sim \mathcal{N}\left(0,K(X_{\ast },X_{\ast })\right)$$从中Sample出三个函数如下图(a)的红、绿、蓝（这些测试点只是为了画出函数）

2. 假设有了观察得到的训练样本点 { x i , f ( x i ) } i = 1 n \{x_i,f(x_i)\}_{i=1}^n {xi​,f(xi​)}i=1n​，拼成训练矩阵$X\in {\mathbb{R}}^{1\times n}$，其输出拼成$\mathbf{f}$，而测试矩阵$X_{\ast }\in {\mathbb{R}}^{1\times n_{\ast }}$，输出为${\mathbf{f}}_{\ast }$，于是它们的joint prior为：
   $$\begin{gathered} \left[\begin{array}{l}\mathbf{f}\\ {\mathbf{f}}_{\ast }\end{array}\right]\sim \mathcal{N}\left(0,\left[\begin{array}{ll}K(X,X)& K\left(X,X_{\ast }\right)\\ K\left(X_{\ast },X\right)& K\left(X_{\ast },X_{\ast }\right)\end{array}\right]\right) \\ \text{w}hereK(X,X_{\ast })\in {\mathbb{R}}^{n\times n_{\ast }},K(X_{\ast },X)\in \mathbb{R}{\in }^{n_{\ast }\times n} \end{gathered}$$

3. 1.2节中的(3)式即Bayesian Linear Regression的推断公式：$$\begin{array}{rl}{f}_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\sim \mathcal{N}(& {\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi {\left(K+{\sigma }_n^{2}I\right)}^{-1}\mathbf{y},\\ & {\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast }-{\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi {\left(K+{\sigma }_n^{2}I\right)}^{-1}{\Phi }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast })\end{array}\text{\hspace{1em}(3)}$$

   此处的${\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi$实际上就是Kernel（嗯嗯！不懂可看文末总结）

   可以得到观察到训练样本后，测试的输出函数分布如图(b)所示:
   $$\begin{array}{rl}{\mathbf{f}}_{\ast }\mid X_{\ast },X,\mathbf{f}\sim \mathcal{N}(& K\left(X_{\ast },X\right)K(X,X{)}^{-1}\mathbf{f}\\ & K\left(X_{\ast },X_{\ast }\right)-K\left(X_{\ast },X\right)K(X,X{)}^{-1}K\left(X,X_{\ast }\right))\end{array}$$
   

• 从函数的prior distribution中Sample出三个函数如图a（蓝、绿、红），这些函数特性由covariance function指定
• 然后根据标有"+"号的样本训练后，这三个函数发生了相应的变化如图b

所以总结一下：

1. 不再像多元高斯分布那样，每一维度对应数据的一个属性，高斯过程把一个具体的数据点（样本）$\vec{\mathbf{x}}$当作一个随机变量，高斯过程回归的目的是预测函数$f(\vec{\mathbf{x}})$的分布
2. 根据高斯过程的定义，有限个随机变量 { x ⃗ 1 , . . . , x ⃗ k } \{\vec x_1,...,\vec x_k\} {x 1​,...,x k​}之间服从联合多元的高斯分布$\mathcal{N}({\mu }_k,{\Sigma }_k)$，而它们的函数值$f({\vec{x}}_i),i\in k$也服从对应的正态分布$\mathcal{N}({\mu }_i,{\Sigma }_i)$
3. 定义高斯过程的prior distribution中covariance function是最为关键的，从输出函数的prior中sample一些函数，观察到observation后根据贝叶斯推断的公式如(3)式得到sampled function的posterior

2.3 更为一般的GPR：

在2.2节中并没有考虑导致测量误差、观测误差的噪声变量，因此更为General的GPR如下概率图所示：

核心的一点是：我们看到的数据是应该带有噪声的
不考虑噪声的数据集在2.2节中表示为 { x i , f ( x i ) ≜ f i } \{\mathbf x_i,f(x_i)\triangleq f_i\} {xi​,f(xi​)≜fi​}，那么考虑噪声后：
$$y=f(\mathbf{x})+ϵ,ϵ\sim \mathcal{N}(0,{\sigma }_n^2)$$

其中${\sigma }_n$表示取决于样本本身的噪声，与其它样本独立，所以带噪声的数据集表示为： { x i , y i ≜ f i + ϵ } \{\mathbf x_i,y_i\triangleq f_i+\epsilon\} {xi​,yi​≜fi​+ϵ}

• 考虑噪声的covariance function表示为：
  $$\mathrm{cov}(\mathbf{y})=K(X,X)+{\sigma }_n^{2}Ior\mathrm{cov}(y_p,y_q)=k({\mathbf{x}}_p,{\mathbf{x}}_q)+{\sigma }_n^2{\delta }_{pq},{\delta }_{pq}\sim \mathcal{N}(0,1)$$

• prior distribution为:
  $$\left[\begin{array}{c}\mathbf{y}\\ {\mathbf{f}}_{\ast }\end{array}\right]\sim \mathcal{N}\left(0,\left[\begin{array}{cc}K(X,X)+{\sigma }_n^{2}I& K\left(X,X_{\ast }\right)\\ K\left(X_{\ast },X\right)& K\left(X_{\ast },X_{\ast }\right)\end{array}\right]\right)$$

• posterior distribution为：
  $$\begin{array}{rl}{\mathbf{f}}_{\ast }\mid X,\mathbf{y},X_{\ast }& \sim \mathcal{N}\left({\overline{\mathbf{f}}}_{\ast },\mathrm{cov}\left({\mathbf{f}}_{\ast }\right)\right),\text{ where }\\ {\overline{\mathbf{f}}}_{\ast }& \triangleq \mathbb{E}\left[{\mathbf{f}}_{\ast }\mid X,\mathbf{y},X_{\ast }\right]=K\left(X_{\ast },X\right){\left[K(X,X)+{\sigma }_n^{2}I\right]}^{-1}\mathbf{y}\\ \mathrm{cov}\left({\mathbf{f}}_{\ast }\right)& =K\left(X_{\ast },X_{\ast }\right)-K\left(X_{\ast },X\right){\left[K(X,X)+{\sigma }_n^{2}I\right]}^{-1}K\left(X,X_{\ast }\right)\end{array}$$

仔细看看mean以及variance的项

1. mean取决于observed targets或观测或label即$\mathbf{y}$、测试点与观测点在高维空间的距离$K(X_{\ast },X)$
2. covariance function与观测值$\mathbf{y}$无关！只与输入有关！方差一般衡量不确定性，所以只与测试点和观测点的距离相关，与直觉相符！
3. covariance function第一项$K(X_{\ast },X_{\ast })$只是单纯的prior covariances没有任何信息，然后减去了一项观测到$X$的信息变化

三、高斯过程总结

3.1 两个角度的联系

用了大量篇幅描述了Weight-Space View以及Function Space View，现在用来点明它们之间的联系，说明它们本质上是同一个东西。
在Weight-Space View中的贝叶斯推断形式(posterior)为：

$$\begin{gathered} \begin{array}{rl}{f}_{\ast }\mid {\mathbf{x}}_{\ast },X,\mathbf{y}\sim \mathcal{N}(& {\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi {\left(K+{\sigma }_n^{2}I\right)}^{-1}\mathbf{y},\\ & {\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast }-{\mathbf{\varphi }}_{\ast }^{\top }{\Sigma }_p\Phi {\left(K+{\sigma }_n^{2}I\right)}^{-1}{\Phi }^{\top }{\Sigma }_p{\mathbf{\varphi }}_{\ast })\end{array} \\ \text{where}K={\Phi }^{\top }{\Sigma }_p\Phi \text{\hspace{1em}(3)} \end{gathered}$$

在Function-Space View中的函数分布推断形式(posterior)为：
$$\begin{array}{rl}{\mathbf{f}}_{\ast }\mid X,\mathbf{y},X_{\ast }& \sim \mathcal{N}\left({\overline{\mathbf{f}}}_{\ast },\mathrm{cov}\left({\mathbf{f}}_{\ast }\right)\right),\text{ where }\\ {\overline{\mathbf{f}}}_{\ast }& \triangleq \mathbb{E}\left[{\mathbf{f}}_{\ast }\mid X,\mathbf{y},X_{\ast }\right]=K\left(X_{\ast },X\right){\left[K(X,X)+{\sigma }_n^{2}I\right]}^{-1}\mathbf{y}\\ \mathrm{cov}\left({\mathbf{f}}_{\ast }\right)& =K\left(X_{\ast },X_{\ast }\right)-K\left(X_{\ast },X\right){\left[K(X,X)+{\sigma }_n^{2}I\right]}^{-1}K\left(X,X_{\ast }\right)\end{array}$$

所以选择Kernel Function $K(C,D)=\Phi (C{)}^{\top }{\Sigma }_p\Phi (D)$，其中$C,D$为样本矩阵或测试矩阵$Xor{X}_{\ast }$，element-wise地看为$k({\mathbf{x}}_p,{\mathbf{x}}_q)=\varphi ({\mathbf{x}}_p{)}^{\top }{\Sigma }_p\varphi ({\mathbf{x}}_q)$

所以高斯过程

1. 从Weight-space view角度看，就是将输入$\mathbf{x}$进行高维映射$\varphi (\mathbf{x})$，需要确定这个$\varphi$，然后这些高维特征与参数是线性关系，建模输出$f(\mathbf{x})$，经过噪声变量$ϵ\sim \mathcal{N}(0,{\sigma }_n)$得到我们的观测或label即$y$，其中参数是一个分布，而不是一个值。
2. 从Function-Space View角度看，是通过Kernel Function设定了高维空间的内积关系$<\varphi (\mathbf{x}),\varphi ({\mathbf{x}}^{\prime })>$，隐式地确定了这个映射$\varphi$，然后建模的是关于输出$f(x)$的分布，实际上也就是参数$w$的分布

3.2 两种View的量化

之前用的都是$K(X,X_{\ast })$这种矩阵描述，训练样本矩阵与测试样本矩阵的距离，有点懵。此处通过一个测试样本点$x_{\ast }$来量化两种View，训练样本点为$x_1,...,x_n$，对应的观测为$y_1,...,y_n$，因此mean和covariance为：

$$\begin{array}{rl}& \overline{f_{\ast }}=\left[\begin{array}{c}k(x_{\ast },x_1)\\ k(x_{\ast },x_2)\\ ⋮\\ k(x_{\ast },x_n)\end{array}\right](K(X,X)+{\sigma }_n^{2}I{)}^{-1}\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]\triangleq {\mathbf{k}}_{\ast }^{\top }(K+{\sigma }_n^{2}I{)}^{-1}\mathbf{y}\\ & \mathbb{V}[f_{\ast }]=k(x_{\ast },x_{\ast })-{\mathbf{k}}_{\ast }^{\top }(K+{\sigma }_n^{2}I{)}^{-1}{\mathbf{k}}_{\ast }\end{array}$$

所以这个测试点$x_{\ast }$的mean输出$\overline{f_{\ast }}$，可理解为训练样本观测值$\mathbf{y}$的线性组合；从kernel的角度看，是n个kernel的线性组合即$\overline{f_{\ast }}={\sum }_i^{n}k(x_{\ast },x_i){\alpha }_i,={\mathbf{k}}_{\ast }^{\top }\alpha$，其中$\alpha =(K+{\sigma }_n^{2}I{)}^{-1}\mathbf{y}$。然后别忘了这个输出值是个后验分布：
$${\mathbf{f}}_{\ast }|X,\mathbf{y},x_{\ast }\sim \mathcal{N}(\overline{f},\mathbb{V}[f_{\ast }])$$

如果还记得贝叶斯：
$$\text{posterior}=\frac{\text{likelihood}\times \text{prior}}{\text{marginal likelihood}}⟺p(\mathbf{w}\mid \mathbf{y},X)=\frac{p(\mathbf{y}|X,\mathbf{w})p(\mathbf{w})}{p(\mathbf{y}\mid X)}$$

有时可以利用marginal likelihood即$p(\mathbf{y}\mid X)$，计算这个的时候有$X\to \mathbf{f}\triangleq f(X)\to \mathbf{y}$，即样本经过模型再经过噪声得到观测的过程。

$$\begin{array}{r}p(\mathbf{y}\mid X)=\int p(\mathbf{y}|\mathbf{f},X)p(\mathbf{f}|X)d\mathbf{f}\end{array}$$

其中$\mathbf{f}|X\sim \mathcal{N}(0,K),\mathbf{y}|\mathbf{f}\sim \mathcal{N}(\mathbf{f},{\sigma }_n^{2}I)$，因此$\mathbf{y}\sim \mathcal{N}(0,K+{\sigma }_n^{2}I)$，所以有:
$$\log p(\mathbf{y}\mid X)=-\frac{1}{2}{\mathbf{y}}^{\top }{\left(K+{\sigma }_n^{2}I\right)}^{-1}\mathbf{y}-\frac{1}{2}\log \left|K+{\sigma }_n^{2}I\right|-\frac{n}{2}\log 2\pi$$

3.3 算法过程

选定高斯过程的covariance function以及噪声变量的noise level，这里用了cholesky分解对矩阵求逆。这个算法是得到了一个测试样本点$x_{\ast }$在样本矩阵$X$以及观测值$\mathbf{y}$下mean以及variance的输出，那一个测试矩阵就同理了。

3.4 Kernel

Kernel有很多种如Linear Kernel，Periodic Kernel，上述介绍的一维度的Squared-exponential covariance function即RBF，公式如下，就有length-scale即$\ell$，signal variance即${\sigma }_f^2$，noise variance即${\sigma }_n^2$三个超参数的调整。

$$\begin{gathered} k_y\left(x_p,x_q\right)={\sigma }_f^2\exp \left(-\frac{1}{2{\ell }^2}{\left(x_p-x_q\right)}^2\right)+{\sigma }_n^2{\delta }_{pq} \\ \text{where}{\delta }_{pq}\sim \mathcal{N}(0,1) \end{gathered}$$

3.5 Decision Rule

一个测试样本$x_{\ast }$，既然通过了3.3节得到了其mean以及variance的函数分布，那选哪个值？就是这里的Decision Rule，一般会选均值，但更为General的是确定一个loss function来进行调整。按照上述posterior得到一个$y_{\ast }$的分布，从该分布中根据Decision Rule选择一个$y_{guess}$，然后根据loss function来从分布中选择最合适的值$y_{optimal}$，降低分布的随机定。

$$\begin{array}{r}{y}_{optimal}|x_{\ast }=\underset{y_{guess}}{\mathrm{argmin}}\int p(y_{\ast }|x_{\ast },X,\mathbf{y})\mathcal{L}(y_{\ast },y_{guess})d{y}_{\ast }\end{array}$$

认真看看这个$y_{optimal}$的选取，就理解了，这个loss实际上就是在衡量不确定性造成的损失，然后从中选出使损失最小的$y_{optimal}$或者说做出最优的决定。

拓展

这个只是最原始的高斯回归过程，其中有很多需要选择的超参数以及改进的地方。
比如算法速度的优化，kernel函数的选取以及影响。具体可参考

Gaussian Process for Machine Learning
visual gaussian process