转自：http://www.datalearner.com/blog/1051459170229238

　　关于高斯过程，其实网上已经有很多中文博客的介绍了。但是很多中文博客排版实在是太难看了，而且很多内容介绍也不太全面，搞得有点云里雾里的。因此，我想自己发表一个相关的内容，大多数内容来自于英文维基百科和几篇文章。

来源1：https://en.wikipedia.org/wiki/Gaussian_process
来源2：Gaussian Processes in Machine Learning
来源3：Gaussian Processes for Regression and Classification: A Quick Introduction

　　我们还是先来个引言吧。随机过程是硕士研究生的课程了，很不幸，我就挂了。然而现在还是无可避免继续接触，只能说命里有时躲不掉啊  。随机过程英文是stochastic process或者random process。在概率和统计中，随机过程是一个统计过程，它可以理解成一组随机变量的变化。简单说就是，有很多点，每个点都是一个随机变量，很多个随机变量点有次序的在一起就成了随机过程。所以可以说随机过程是一个关于时间t（连续区间可以是时间也可以是空间之类的）分布。每个t都有个对应的随机变量。当t是离散的时候，它就是一组随机变量的序列，比如markov chain。那么高斯过程就可以理解成每个t的值都是来自于一个高斯分布。高斯过程很重要，原因是它已经被研究的很成熟了，如果我们把问题用高斯过程建模，那么很多相关知识我们就已经知道了。高斯过程的应用有很多，可以作为贝叶斯推断的先验，也可以用来做回归和分类。
　　我们用如下的R语言代码生成了一个二维高斯分布。

 
 
x<-y<-seq(-4,4,length=40)
f<-function(x,y){(exp(-0.5*x^2-0.5*y^2))/(2*pi)}
z<-outer(x,y,f)
persp(x,y,z,theta=45,phi=35)

---

1、高斯过程

　　定义：高斯过程是指随机变量的一个集合，其中任意有限个样本的线性组合都有一个联合高斯分布。

　　一个高斯过程是由均值函数$m(x)$和协方差函数k(x,x^{*})k(x,x​∗​​)确定的。它可理解成高斯分布的一个生成过程（或者叫泛化？generalization）。高斯分布的均值和协方差是向量和矩阵（意思就是多维高斯分布的均值和方差是确定的值，比如均值是(0.1,0.2,0.5)。），而高斯过程的均值和方差则分别是均值函数和协方差矩阵函数。我们有如下公式（1）：

$f\sim gp(m,k)$

　　它的含义是函数$f$是一个高斯过程，其均值函数为$m$，协方差函数为$k$。
　　我们这里将简单说一下从分布到随机过程的生成。其实就是说明一下怎么把高斯过程理解成一个分布。来自高斯分布的一个向量中的单个随机变量的索引是它们在向量中的位置。对于高斯过程来说，随机函数$f(x)$的参数$x$就是索引集合，每个输入的$x$都与随机变量$f(x)$关联，也就是这个函数在这个位置上的值。为了方便表示，我们举个例子，假设$x$的取值是自然数（也就是$1,2,..,n$），同时这些值也就是索引。考虑如下的高斯过程：

$f\sim gp(m,k)$

　　其中，m(x)=\frac{1}{4}x^{2}m(x)=​4​​1​​x​2​​，k(x,x^{*})=exp(-\frac{1}{2}(x-x^{*})^{2})k(x,x​∗​​)=exp(−​2​​1​​(x−x​∗​​)​2​​)，(2)

　　为了理解这个过程，我们可以从函数$f$中抽取样本。为了只处理有限情况，我们要求$f$是有限的。我们如何抽取样本呢？给定$x$的值我们可以利用公式（2）来估算均值向量和协方差矩阵（下面是公式3）。

\mu_{i} = m(x_{i}) = \frac{1}{4}x^{2}, i=1,...,nμ​i​​=m(x​i​​)=​4​​1​​x​2​​,i=1,...,n

\sum_{ij} = k(x_{i},x_{j}) = exp(-\frac{1}{2}(x_{i}-x_{j})^{2}),i,j=1,...,n.∑​ij​​=k(x​i​​,x​j​​)=exp(−​2​​1​​(x​i​​−x​j​​)​2​​),i,j=1,...,n.

　　为了区别随机过程和分布，我们使用$m$和$k$表示过程，$\mu$和$\sum$表示分布。我们现在可以从这个分布中产生一个向量（公式4）。

$f\sim N(\mu ,\sum )$

　　我们可以使用matlab画出$f$的图像（如图1所示）：

 
 
xs = (-5:0.2:5)’; ns = size(xs,1); keps = 1e-9;
m = inline(’0.25*x.^2’);
K = inline(’exp(-0.5*(repmat(p’’,size(q))-repmat(q,size(p’’))).^2)’);
fs = m(xs) + chol(K(xs,xs)+keps*eye(ns))’*randn(ns,1);
plot(xs,fs,’.’)

　　在这个例子中，$m$和$k$是均值和协方差，chol是用来计算矩阵的Cholesky分解的（这个分解我也不知道是啥，(◎﹏◎)）。

2、后验高斯过程

　　在前面的例子中，我们说明了如何使用高斯过程定义函数上的分布。高斯过程可以用来作为贝叶斯推断的先验。先验不需要依赖于任何训练数据，但却可以确定函数的一些属性。举个例子来说，图1的函数是平滑的，且与二次方程的图像很接近。这部分内容的目标就是要在包含训练数据的情况下如何推导出简单的规则来更新这个先验（这里可以解释更新先验的概念，其实就是从数据中引入信息到先验中，然后把我们的先验信念转化成后验信念。开始假设先验的时候都是独立于数据，但与实际上可能是有参数上的差距【因为模型如果有差距是测不出来的】，如果我们计算发现后验概率与先验概率差别较大，那么可能需要我们调整先验参数使其更加符合实际，所以这里更新先验可以理解成计算后验概率，可参见http://stats.stackexchange.com/questions/166321/what-exactly-does-it-mean-to-and-why-must-one-update-prior 的回答：In plain english, update a prior in bayesian inference means that you start with some guesses about the probability of an event occuring (prior probability), then you observe what happens (likelihood), and depending on what happened you update your initial guess. Once updated, your prior probability is called posterior probability.）。下一节我们将借助于数据说明这个先验的一些属性。这里的意思就是，假设先验是高斯过程，我们如何利用训练数据来更新这个先验的分布。
　　后验函数可以用来预测新的数据。假设f是通过训练数据已知的函数值，f _{*}​∗​​ 是测试集输入X_{*}​∗​​的函数值。现在我们可以把我们感兴趣的变量写成联合分布：

　　我们说明一下符号：\mu=m(x_{i}),i=1,...,nμ=m(x​i​​),i=1,...,n是训练集的均值。测试集均值为\mu_{*}μ​∗​​。同样的，训练集的协方差我们使用\sum_{*}∑​∗​​，测试集的协方差为\sum_{**}∑​∗∗​​。假设训练集的f已知，我们感兴趣的是f_{*}​∗​​的条件概率（公式6）：

f_{*}|f \sim N(\mu_{*}+\sum^{T}_{*}\sum^{-1}(f-\mu),\sum_{**}-\sum^{T}_{*}\sum^{-1}\sum_{*})f​∗​​∣f∼N(μ​∗​​+∑​∗​T​​∑​−1​​(f−μ),∑​∗∗​​−∑​∗​T​​∑​−1​​∑​∗​​)

　　这是特定测试集下的后验分布。很容易验证相应的后验过程是（公式7）：

f|D \sim gp(m_{D},k_{D})f∣D∼gp(m​D​​,k​D​​)

m_{D}(x) = m(x)+\sum(X,x)^{T}\sum^{-1}(f-m)m​D​​(x)=m(x)+∑(X,x)​T​​∑​−1​​(f−m)

k_{D}(x,x') = k(x,x')-\sum(X,x)^{T}\sum^{-1}\sum(X,x')k​D​​(x,x​′​​)=k(x,x​′​​)−∑(X,x)​T​​∑​−1​​∑(X,x​′​​)

　　其中，$\sum (X,x)$是训练集与$x$的协方差向量。这些就是高斯过程预测的核心公式。我们来看一下后验均值和方差。注意到后验方差k_{D}(x,x)k​D​​(x,x)等于先验$k(x,x)$减去一个正值，改正直依赖于训练数据的输入。因此后验方差总是小于先验方差，原因是数据给了我们额外的信息（信息越多，方差应该越小）。
　　同时，我们还是有一个问题需要解决，即训练集的噪音问题。在回归问题中有很多观测量是有噪音的。最常见的假设是它们是独立同分布的。输出具有高斯噪音。在高斯过程模型中，这样的噪音很容易考虑。即每一个$f(x)$都有个额外的协变量，其大小等于噪音的方差（以下为公式8）：

y(x)=f(x)+\varepsilon,\varepsilon\sim N(0,\sigma_{n}^{2})y(x)=f(x)+ε,ε∼N(0,σ​n​2​​)

f \sim gp(m,k),y\sim gp(m,k+\sigma^{2}_{n}\delta_{ii'})f∼gp(m,k),y∼gp(m,k+σ​n​2​​δ​ii​′​​​​)

　　其中，\delta_{ii^*}=1δ​ii​∗​​​​=1当且仅当i=i^*i=i​∗​​是Kronecker delta（这个还未搞清楚啥东西）。注意到Kronecker delta的索引是位置，不是输入量$x$。而不同的输入之间的噪音是相互独立的。因此，噪音过程的协方差函数是信号协方差与噪音协方差的和。至此，我们可以把后验协方差函数画出来，如图2所示。但到这里，我们依然还有些问题没有解决：1）如何找出先验中均值和协方差的函数？2）如何估计噪音的水平？这将在下面一章节中讲述。

３、高斯过程的训练

　　在机器学习中，先验信息通常并不知道。为了使GP具有使用价值，我们必须能够根据数据在不同的均值和协方差之间做出选择。
　　通常我们选择层次先验。即均值和协方差来自于超参数决定的分布（函数）中。这也是贝叶斯分析的基本原理了。比如，公式2中，我们可以把均值函数和协方差函数定义成如下参数决定的函数（公式9）：

$f\sim gp(m,k)$

m(x)=ax^{2}+bx+c, and, k(x,x')=\sigma^{2}_{y}exp(-\frac{(x-x')^{2}}{2l^{2}})+\sigma^{2}_{n}\delta_{ii'}m(x)=ax​2​​+bx+c,and,k(x,x​′​​)=σ​y​2​​exp(−​2l​2​​​​(x−x​′​​)​2​​​​)+σ​n​2​​δ​ii​′​​​​

　　这里的超参数是\theta={a,b,c,\sigma_{y},\sigma_{n},l}θ=a,b,c,σ​y​​,σ​n​​,l。分层可以帮助我们把模糊的先验信息以简单的方式表达。举个例子来说，假设我们认为函数是二阶多项式，但我们不知道它的参数。实际上，二阶多项式与数据之间的差别在于光滑曲线加上噪音函数，但我们不需要指定特征长度范围$l$，也不需要知道二者的量级。我们可以根据数据来推断所有的超参数。
　　为了计算超参数，我们需要计算在超参数条件下数据的概率。幸运的是，这很简单，假设数据分布式高斯分布（公式10）：

L=logp(y|x,\theta)=-\frac{1}{2}log|\sum|-\frac{1}{2}(y-\mu)^{T}\sum^{-1}(y-\mu)-\frac{n}{2}log(2\pi)L=logp(y∣x,θ)=−​2​​1​​log∣∑∣−​2​​1​​(y−μ)​T​​∑​−1​​(y−μ)−​2​​n​​log(2π)

　　我们把这个量乘坐是边缘似然（marginal likelihood）的log函数。我们使用“边缘”来强调我们处理的是非参数模型。现在我们可以根据边缘似然的偏导计算超参数的值（公式11）。

\frac{\partial L}{\partial \theta_{m}} = -(y-\mu)^{T}\sum^{-1}\frac{\partial m}{\partial\theta_{m}}​∂θ​m​​​​∂L​​=−(y−μ)​T​​∑​−1​​​∂θ​m​​​​∂m​​

\frac{\partial L}{\partial \theta_{k}} = \frac{1}{2}trace(\sum^{-1}\frac{\partial\sum}{\partial \theta})+\frac{1}{2}(y-\mu)^{T}\frac{\partial \sum}{\partial \theta_{k}}(y-\mu)​∂θ​k​​​​∂L​​=​2​​1​​trace(∑​−1​​​∂θ​​∂∑​​)+​2​​1​​(y−μ)​T​​​∂θ​k​​​​∂∑​​(y−μ)

　　其中\theta_{m}θ​m​​和\theta_{k}θ​k​​分别是超参数的均值和方差函数。这个公式的求解还有点不太懂，但是该文作者说Eq. (11) can conveniently be used in conjunction with a numerical optimization routine such as conjugate gradients to find good hyperparameter settings.还是数学太差没办法啊。（注意，这个函数不是凸函数，不能简单求极值）。
　　由于高斯过程是非参数模型，其边缘似然与正常的参数模型不太一样。注意到，模型适配数据非常简单，只要我们把噪音水平\sigma^{2}_{n}σ​n​2​​设置成0，模型就会输出均值预测函数，那么它与训练集就完全吻合了。但这并不是正确的优化边缘似然的方法。实际上，公式10有三个部分：第一个是复杂性惩罚量（complexity penality term），用来测度并降低模型复杂性的；第二个是负的多项式函数，作用是适配数据测量的（是唯一一个依赖于训练集的）；第三个是Log的正则项，独立于数据之外，作者说是not interesting。图3是用极大边缘似然得到的预测模型。
　　注意到，惩罚项和数据适配在GP模型中是自动的。我们不需要额外的方法来确定权重参数。这在实际中非常重要，因为它简化了训练。图4说明了这种均衡是如何出来的。

---

4、如何使用高斯过程解决回归

首先来看另一种高斯过程的解释
　　高斯过程将多元高斯分布扩展到无限维度。任意数据集中n个观测量，y={y_{1},..,y_{n}}y=y​1​​,..,y​n​​都可以理解成从某个多元高斯分布中抽取的一个样本。通常，我们都将高斯过程的均值设置为0，协方差最常用的选择是“平方指数（squared exponential）”：

k(x,x')=\sigma^{2}_{f}\exp[\frac{-(x-x')^{2}}{2l^2}]k(x,x​′​​)=σ​f​2​​exp[​2l​2​​​​−(x−x​′​​)​2​​​​]

　　由于噪音的存在，我们可以把模型建立成如下的形式：

y=f(x)+N(0,\sigma^{2}_{n})y=f(x)+N(0,σ​n​2​​)

k(x,x')=\sigma^{2}_{f}\exp[\frac{-(x-x')^{2}}{2l^2}]+\sigma^{2}_{n}\delta(x,x')k(x,x​′​​)=σ​f​2​​exp[​2l​2​​​​−(x−x​′​​)​2​​​​]+σ​n​2​​δ(x,x​′​​)

　　其中，\delta(x,x^{*})δ(x,x​∗​​)是KroneckerDelta函数（这里其实就是我们常说的指示函数的意思，当元素二者相等时取1，否则取0）。为了获得高斯过程回归GPR，我们需要计算上述协方差函数。这里的x_{*}x​∗​​是指当前位置，也就是新观测量的自变量。因为我们要预测在当前训练集下新自变量的输出值，所以计算的时候除了需要训练集的协方差矩阵外，我们还需要新的自变量与训练集自变量之间的协方差，计算出下面三个矩阵就可以方便我们后面的计算了。

　　由于高斯过程的可以把数据当作从多元高斯分布的的一个样本，所以我们有：

　　与前面相同，我们可以计算出边缘似然为公式（7）

y_{*}|y \sim N(K_{*}K^{-1}y,K_{**}-K_{*}K^{-1}K^{T}_{*})y​∗​​∣y∼N(K​∗​​K​−1​​y,K​∗∗​​−K​∗​​K​−1​​K​∗​T​​)

　　那么y_*y​∗​​的均值和方差估计公式为：

\bar{y}_{*}=K_{*}K^{-1}y​y​¯​​​∗​​=K​∗​​K​−1​​y

var(y_{*})=K_{**}-K_{*}K^{-1}K^{1}_{*}var(y​∗​​)=K​∗∗​​−K​∗​​K​−1​​K​∗​1​​

　　让我们举个例子，如下图所示：

　　我们有6个观测值y，分别位于：

$x=[-1.50,-1.00,-0.75,-0.40,-0.25,-0.00]$

　　我们从误差线得知\sigma_{n}=0.3σ​n​​=0.3，当我们审慎而明智的选择了\sigma_{f}σ​f​​和$l$的值的时候（作者用词真是厉害，审慎而明智的是with judicious，怎么选择后面会说）我们就可以计算协方差矩阵了。

　　我们知道了K_{**}=1.70K​∗∗​​=1.70，也就知道了

K_{*}=[0.38,0.70,1.03,1.35,1.46,1.58]K​∗​​=[0.38,0.70,1.03,1.35,1.46,1.58]

　　进而我们知道了\bar{y}_{*}=0.95​y​¯​​​∗​​=0.95，var(y_{*})=0.21var(y​∗​​)=0.21。如下图所示：