研究

什么是高斯过程

在讨论高斯过程(GP)之前,我们先从最简单的高斯分布开始。高斯分布是一种概率分布,相信学过概率论的人都了解什么是概率分布。

从高斯分布开始

在讨论高斯过程(GP)之前,我们先从最简单的高斯分布开始。高斯分布是一种概率分布,相信学过概率论的人都了解什么是概率分布。

对于简单的一元高斯分布,其定义如下:

p(x)=1σ2πexp((xμ)22σ2)p(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp (-\frac{(x-\mu)^2}{2\sigma^2})

在给定均值 μ\mu 和方差 σ\sigma 之后,我们就可以确定一个高斯分布,并计算出事件 xx 所对应的概率。

如果我们希望计算两个事件 x1,x2x_1,x_2 的高斯分布呢?更一般的,如何计算多个事件的高斯分布呢?这就需要将一元高斯分布推广到多元,假设多元高斯分布的各个维度之间相互独立,那么我们可以推导出以下公式:

p(x1,x2,,xn)=i=1np(xi)=1(2π)n2σ1σ2σnexp(12[(x1μ1)2σ12+(x2μ2)2σ22++(xnμn)2σn2])\begin{aligned} p(x_1,x_2,\dots,x_n)&=\prod_{i=1}^{n}p(x_i)\\ &= \frac{1}{(2\pi)^{\frac{n}{2}}\sigma_1\sigma_2\dots\sigma_n}\exp\left(-\frac{1}{2}\left[\frac{(x_1-\mu_1)^2}{\sigma_1^2}+\frac{(x_2-\mu_2)^2}{\sigma_2^2}+\dots+\frac{(x_n-\mu_n)^2}{\sigma_n^2}\right]\right) \end{aligned}

这是一个比较简单的公式推导,通过概率论的乘法法则定义了多个事件下的多元高斯分布,即多个一元高斯分布的连乘。

对于连乘连加我们应当自然而然的想到线性代数,因此接下来我们通过线性代数简化一下上述的表示。 xμ=[x1μ1,x2μ2,,xnμn]TK=[σ12000σ2200000]\begin{aligned} \mathbf{x}-\mathbf{\mu}&=\left[x_1-\mu_1,x_2-\mu_2,\dots,x_n-\mu_n\right]^T \\ \mathbf{K}&= \begin{bmatrix} \sigma_1^2 & 0 & \dots & 0\\ 0 & \sigma_2^2 & \dots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 0 & 0 & \cdots & 0 \end{bmatrix} \end{aligned}

为什么要这样做?我们需要注意,左侧分式的分母是连乘,这自然是一个对角阵。而右侧分式则是元素乘法求和,这自然是一个向量和对角阵的值。

没错,我也不知道怎么冒出来的,可能数学就是经验吧。

由此我们可以进行如下简化: σ1σ2σn=K12(x1μ1)2σ12+(x2μ2)2σ22++(xnμn)2σn2=(xμ)TK1(xμ)\begin{aligned} \sigma_1\sigma_2\cdots\sigma_n&=\lvert\mathbf{K}\rvert^{\frac{1}{2}}\\ \frac{(x_1-\mu_1)^2}{\sigma_1^2}+\frac{(x_2-\mu_2)^2}{\sigma_2^2}+\dots+\frac{(x_n-\mu_n)^2}{\sigma_n^2}&=(\mathbf{x}-\mathbf{\mu})^T\mathbf{K}^{-1}(\mathbf{x}-\mathbf{\mu}) \end{aligned}

你问我怎么写出来,说真的我也不知道,因为这不是我写的。

代入多元高斯分布公式得:

p(x)=(2π)n2K12exp(12(xμ)TK1(xμ))p(x)=(2\pi)^{-\frac{n}{2}}\lvert\mathbf{K}^{-\frac{1}{2}}\rvert\exp\left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu})^T\mathbf{K}^{-1}(\mathbf{x}-\mathbf{\mu})\right)

其中 μ\mathbf{\mu} 为均值向量,K\mathbf{K} 为协方差矩阵(在各维度独立下为对角阵)。实际上这个简化后的公式就是线代版的高斯分布公式。因此也可以简写为:

xN(μ,K)x\sim \mathcal{N} \left( \mathbf{\mu}, \mathbf{K}\right)

如果各个事件变量不独立,则 K\mathbf{K} 不为对角阵,形式不变。

什么是高斯过程

前面我们定义了多元高斯分布,进一步说,实际上我们定义了对于多个随机变量的高斯分布。而所谓高斯过程,则是关于随机过程的高斯分布。

那么什么是随机过程呢?随机过程就是发生在一个连续域中的随机变量的集合。这听起来不明所以,我们举个例子。我们定义每间隔一个小时发生一个事件,这个事件发生时遵循一个高斯分布。

我们这里假设遵循高斯分布,直接上就是直接定义了一个高斯过程。

没错,这个不断发生的事件遵循什么分布就是什么过程。

例如我假设每间隔一个小时,你正在看的这个屏幕会有有可能爆炸, 而发生爆炸的概率遵循一个高斯分布(注意,前一个小时爆炸的概率所遵循的高斯分布和后一个小时所遵循的高斯分布不需要一致)。因此当时间不断流逝,我们相当于在每个小时的高斯分布上进行了一次取样,我们把所有的取样连接起来,构成了一条曲线。

现在假设我们不幸遇到了爆炸,我们有机会从头重新开始经历这个过程(不爆炸也可以)。即我们重新在时间轴上,对每个时刻重新进行一次采样,我们又得到了一条曲线。不断重复这个过程,我们就可以得到很多曲线。

现在我们从函数的角度,认为每一条曲线代表了一个函数,因此我们实际上是在不断取样得到新的函数,整个过程描述了我们所采样得到的函数的分布状态,因此高斯过程又被称为之函数的高斯分布。

现在你也许有点摸不着头脑,到底什么是高斯过程?实际上我们所定义的 每间隔一个小时发生一个事件,这个事件发生时遵循一个高斯分布 就是所谓的高斯过程。即在一个无限连续域(时间)上的服从高斯分布的随机变量(爆炸)的集合。也被称为无限高斯分布。

高斯过程还有一个前提,就是连续域上任意有限个的随机变量的联合分布遵循一个多元高斯分布。

需要注意,高斯过程意味着所有变量的边缘分布为高斯分布,但是所有变量遵循高斯分布不代表整体是一个高斯过程。

我们自然而然可以想到,如果我们采样均值和方差,那么就可以得到一个均值函数 m(t)m(t) 和一个协方差函数 k(t)k(t)。这里的协方差函数假定了每个时刻的随机变量相互独立,通常而言我们会认为不同时刻之间是相互影响的,因此我们将其推广为 k(s,t)k(s,t),表示时刻 s 和时刻 t 的随机变量的协方差函数,又被称为核函数。

协方差函数是一种核函数,不要理解为核函数就是协方差函数。

这里由 k(t)k(t) 推广出的 k(s,t)k(s,t) 实际上描述了一个协方差矩阵,对于不同的 s 和 t 取值,这个函数可以生成一个矩阵用于描述所有时刻之间的方差,而对角线就是某一时刻其本身的方差。

接下来我们用数学定义一下高斯过程,使用 ξt\xi_t 表示时间 t 时刻的随机变量,它服从一个均值为 m(t)m(t) 方差为 k(s,t)k(s,t) 的高斯分布。这里的 s 表示另一个时刻 s。此时我们就确定了一个高斯过程 fGP(m(t),k(s,t))f\sim GP(m(t),k(s,t))。这是一个描述了函数的分布。

这里我们引入了其他时刻对于 t 时刻的影响,是因为在应用高斯过程的时候,通常我们认为不同时刻之间会相互影响。实际上如果不同时刻之间的高斯分布相互独立,那么 k(s,t)k(s,t) 就退化为 k(t)k(t),这虽然符合我们对于一元高斯分布的认知,但是没有任何用处。因为这代表着一个服从高斯分布的噪声,前后没有任何关联。

高斯过程回归

对于一个高斯过程 fGP(m(t),k(s,t))f\sim GP(m(t),k(s,t)),我们可以像做一元高斯回归一样,对高斯过程进行回归,从而确定其均值函数和协方差函数。

实际上我们也可以理解高斯过程回归是一种函数回归方法,但是不同于点估计,例如线性回归等。高斯过程回归是一种贝叶斯估计,高斯回归的结果不会给出一个明确的函数表达式,或者说给出一个非参数化的结果。而且这个结果在给出预测的同时,会给出这个预测所服从的分布,从而我们可以确定预测值的置信度。

所谓参数化,即假定模型有固定数量的参数,例如神经网络,从一开始参数的数量就固定了。而高斯过程的预测依赖于协方差函数,协方差函数生产的协方差矩阵会随着数据量的增多而变大,因此其参数量是随着数据量变化的,称之为非参数化。

假设我们有一些观测值 x_train 和 y_train。现在我们希望使用这些数据进行高斯过程回归,从而可以对新数据进行预测。我们需要做什么呢?

为什么要使用高斯过程回归?我不可以直接用线性回归或者多项式回归或者神经网络吗?

当然可以,你可以选择任意回归方式,只是不同的回归方式其效果也不同。例如我们不可能对一个曲线采用线性回归。

那什么时候我们可以使用高斯过程回归?通常而言,数据量不大(因为高斯过程回归的时间复杂度为 O(n3)O(n^3)),而我们预测不确定性,我们预想中的结果具有一定的平滑性等特点。此时我们就可以选择高斯过程回归。

以下内容生成自 AI

1. 问题设定

假设我们观测到一组输入输出数据对 D={(xi,yi)}i=1n\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n,其中数据由如下噪声模型给出:

yi=f(xi)+εi,εiN(0,σn2)y_i = f(x_i) + \varepsilon_i, \quad \varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)

我们的任务是:在给定观测数据的基础上,对新的输入 xx_* 预测其对应输出 yy_* 的分布。


2. 高斯过程作为分布的先验

我们假设目标函数 f(x)f(x) 是一个高斯过程(Gaussian Process, GP):

f(x)GP(m(x),k(x,x))f(x) \sim \mathcal{GP}(m(x), k(x,x'))
  • m(x)m(x) 为均值函数,常设为0。
  • k(x,x)k(x,x')为核函数,例如 RBF(高斯)核。

含义说明:对任意有限一组输入点 (x1,,xn)(x_1, \ldots, x_n) ,其对应的函数值 (f(x1),,f(xn))(f(x_1),\ldots, f(x_n))联合分布是高斯分布。


3. 写出关联的联合分布

记观测输入列向量 X=[x1,,xn]T\mathbf{X} = [x_1, \ldots, x_n]^T,观测输出 y=[y1,,yn]T\mathbf{y} = [y_1, \ldots, y_n]^T

记待预测点为 xx_*f=[f(x1),,f(xn)]Tf=f(x)\begin{aligned} \mathbf{f} &= [f(x_1), \ldots, f(x_n)]^T\\ f_* &= f(x_*) \end{aligned} 按照高斯过程的定义,有:

[ff]N([mm],[K(X,X)K(X,x)K(x,X)K(x,x)])\begin{bmatrix} \mathbf{f} \\ f_* \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} \mathbf{m} \\ m_* \end{bmatrix}, \begin{bmatrix} K(\mathbf{X},\mathbf{X}) & K(\mathbf{X},x_*) \\ K(x_*,\mathbf{X}) & K(x_*,x_*) \end{bmatrix} \right)

其中

  • m=[m(x1),,m(xn)]T\mathbf{m} = [m(x_1), \ldots, m(x_n)]^T,常为0
  • m=m(x)m_* = m(x_*),常为0
  • K(X,X)K(\mathbf{X},\mathbf{X})n×nn\times n 的核矩阵,
  • K(X,x)=[k(x1,x),,k(xn,x)]TK(\mathbf{X},x_*) = [k(x_1,x_*), \ldots, k(x_n, x_*)]^Tn×1n \times 1 列向量
  • K(x,x)=k(x,x)K(x_*,x_*) = k(x_*,x_*)

解释:这反映了GP对函数空间的协方差建模,仅由核函数控制。


4. 加入观测噪声

我们实际观测到的是 y\mathbf{y},而非真值 f\mathbf{f}

y=f+ε,εN(0,σn2I)\mathbf{y} = \mathbf{f} + \boldsymbol{\varepsilon}, \quad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0}, \sigma_n^2 I)

因此推断时须对训练点加上观测噪声,核矩阵需加一项:

Ky:=K(X,X)+σn2IK_y := K(\mathbf{X},\mathbf{X}) + \sigma_n^2 I

5. 条件高斯分布的推导

(1) 联合分布重写

观测向量 y\mathbf{y} 和目标 ff_*,它们的联合分布是

[yf]N([00],[KyK(X,x)K(x,X)k(x,x)])\begin{bmatrix} \mathbf{y} \\ f_* \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} K_y & K(\mathbf{X},x_*) \\ K(x_*,\mathbf{X}) & k(x_*,x_*) \end{bmatrix} \right)

解释:对观测点加入噪声,所以对角元素加 σn2\sigma_n^2


(2) 利用高斯条件分布公式

设联合分布

[ab]N([μaμb],[ACCTB])\begin{bmatrix} \mathbf{a} \\ \mathbf{b} \end{bmatrix} \sim \mathcal{N}\left( \begin{bmatrix} \boldsymbol\mu_a \\ \boldsymbol\mu_b \end{bmatrix}, \begin{bmatrix} A & C \\ C^T & B \end{bmatrix} \right)

则,ba\mathbf{b}|\mathbf{a} 的条件分布为:

  • E[ba]=μb+CTA1(aμa)\mathbb{E}[\mathbf{b}|\mathbf{a}] = \boldsymbol{\mu}_b + C^T A^{-1} (\mathbf{a} - \boldsymbol{\mu}_a)
  • Cov[ba]=BCTA1C\mathrm{Cov}[\mathbf{b}|\mathbf{a}] = B - C^T A^{-1} C

(3) 套用公式

  • a=y\mathbf{a} = \mathbf{y}, b=f\mathbf{b} = f_*
  • μa=μb=0\boldsymbol\mu_a = \boldsymbol\mu_b = 0
  • A=KyA = K_y, C=K(X,x)C = K(\mathbf{X}, x_*), B=k(x,x)B = k(x_*,x_*)

则:

预测均值:

E[fy]=K(x,X)Ky1y\mathbb{E}[f_* | \mathbf{y}] = K(x_*,\mathbf{X}) K_y^{-1} \mathbf{y}

(常写成 KTKy1yK_*^T K_y^{-1} \mathbf{y},其中 K=K(X,x)K_* = K(\mathbf{X}, x_*)

预测方差:

Var[fy]=k(x,x)K(x,X)Ky1K(X,x)\mathrm{Var}[f_* | \mathbf{y}] = k(x_*, x_*) - K(x_*,\mathbf{X}) K_y^{-1} K(\mathbf{X},x_*)

解释

  • 第一项是测试点在先验下的方差,减去“数据约束”后的修正。
  • 假如输入点 xx_* 离观测点越远,修正项趋向0,预测方差趋于先验
  • Ky1yK_y^{-1} \mathbf{y} 是“数据对预测的权重”,由各点与测试点的核相关性决定

6. 对输出的观测预测

如果感兴趣的是真正的观测 yy_*,要加上噪声:

y=f+ε,εN(0,σn2)y_* = f_* + \varepsilon_*, \quad \varepsilon_* \sim \mathcal{N}(0, \sigma_n^2)

则有

  • E[yy]=E[fy]\mathbb{E}[y_*|\mathbf{y}] = \mathbb{E}[f_*|\mathbf{y}]
  • Var[yy]=Var[fy]+σn2\mathrm{Var}[y_*|\mathbf{y}] = \mathrm{Var}[f_*|\mathbf{y}] + \sigma_n^2

7. 超参数优化:边缘似然函数

核函数(如RBF)与噪声方差有超参数 θ={σf,l,σn}\theta = \{\sigma_f, l, \sigma_n\}

边缘似然(对数)为:

logp(yX,θ)=12yTKy1y12logKyn2log(2π)\log p(\mathbf{y} | \mathbf{X}, \theta) = -\frac{1}{2} \mathbf{y}^T K_y^{-1} \mathbf{y} - \frac{1}{2} \log|K_y| - \frac{n}{2} \log (2\pi)

优化此函数可调优 θ\theta


8. 推导流程小结

  1. 写出高斯过程联合分布,将测试点加入观测点
  2. 用噪声模型调整核矩阵(加 σn2I\sigma_n^2 I
  3. 用高斯条件分布公式推出测试点的后验均值、方差
  4. 对超参数可用边缘似然最大化学习

9. 物理/直观解释

  • 均值公式:每个观测点对预测的“影响力”由核相关度决定(距离远相关小);
  • 方差公式:观测集约束越多预测方差减小,数据多分布密集之处模型更“自信”;
  • 噪声控制:噪声越大核矩阵主对角线越大,表示对数据不太信任,学习更保守。

10. 符号快速总结

观测:y=[y1,,yn]T核矩阵:Ky=K(X,X)+σn2I测试点核向量:K=K(X,x)均值:μ=KTKy1y方差:σ2=k(x,x)KTKy1K\begin{align*} \text{观测:} &\quad \mathbf{y} = [y_1, \ldots, y_n]^T \\ \text{核矩阵:} &\quad K_y = K(\mathbf{X},\mathbf{X}) + \sigma_n^2 I \\ \text{测试点核向量:} &\quad K_* = K(\mathbf{X},x_*) \\ \text{均值:} &\quad \mu_* = K_*^T K_y^{-1} \mathbf{y} \\ \text{方差:} &\quad \sigma^2_* = k(x_*,x_*) - K_*^T K_y^{-1} K_* \end{align*}

References

  1. 高斯过程 Gaussian Processes 原理、可视化及代码实现
  2. 机器学习-白板推导系列(二十)-高斯过程GP(Gaussian Process)