【参考资料】 1.B站：机器学习-白板推导系列(二十)-高斯过程GP(Gaussian Process) 2.：一文搞定（二）—— 贝叶斯线性回归 3.知乎：Gaussian process regression的导出——权重空间视角下的贝叶斯的方法 4.知乎：高斯过程 Gaussian Processes 原理、可视化及代码实现

【文章结构】 1 Introduction 简要介绍 2 Weigh-space view GPR 的权重空间角度 3 Fron wsv to fsv 从权重空间到函数空间是如何过渡的 4 Function-space view GPR 的函数空间角度 5 Example in code 包含对于一个简单 GPR 示例的 python 和 matlab 两种语言的代码

【阅读建议】 对于希望直接应用的读者可以跳过 Section 2, 3 两部分，直接看 Section 1, 4, 5

1 Introduction

1.1 多元高斯分布

由多元高斯分布的定义可知： $$p(\mathbf{x})=\prod _{i=1}^{n}p(x_i)=\frac{1}{(2\pi {)}^{\frac{n}{2}}{\sigma }_1...{\sigma }_n}\exp (-\frac{1}{2}[\frac{(x_1-{\mu }_1{)}^2}{{\sigma }_1^2}...\frac{(x_n-{\mu }_n{)}^2}{{\sigma }_2^2}])\text{\hspace{1em}(1)}$$ 进一步的，令 $\mathbf{x}-\mathbf{\mu }=[x_1-{\mu }_1,...,x_n-{\mu }_n{]}^T,\Sigma =\left[\begin{array}{cccc}{\sigma }_1^2& 0& ...& 0\\ 0& {\sigma }_2^2& ...& 0\\ ...& ...& ...& ....\\ 0& ...& 0& {\sigma }_n^2\end{array}\right]$，从而得到式 $(2)$：

$$p(\mathbf{x})=(2\pi {)}^{-\frac{n}{2}}|\Sigma {|}^{-\frac{1}{2}}\exp (-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^T|\Sigma {|}^{-1}(\mathbf{x}-\mathbf{\mu }))\text{\hspace{1em}(2)}$$

我们将式 $(2)$ 称之为多元高斯分布的向量化表示。注意：留意上式中的 $(\mathbf{x}-\mathbf{\mu }{)}^T|\Sigma {|}^{-1}(\mathbf{x}-\mathbf{\mu })$，它与下文将要介绍的 kernel function 有着类似的形式

1.2 Gaussain Process 高斯过程

高斯过程是定义在连续域上的无限多个高维随机变量所组成的随机过程，可以看做是一个无限维的高斯分布

数学定义

对于随机变量 $\{{\zeta }_t{\}}_{t\in T}$, $T$ 可以是连续的时间或空间

如果对于 $\forall n\in {\mathbb{N}}^{+},t_1,...,t_n\to {\zeta }_1,...,{\zeta }_n$, 都有 $\mathbf{\zeta }=({\zeta }_1,...,{\zeta }_n{)}^T\sim N(\mathbf{\mu },\mathbf{\Sigma })$,

那么 $\{{\zeta }_t{\}}_{t\in T}$ 就是一个高斯过程，记做 $${\zeta }_t\sim GP(m(t),\kappa (s,t))\text{\hspace{1em}(3)}$$

where,

$m(t)=E[{\zeta }_t]$, mean function$\kappa (t,s)=E[({\zeta }_s-m(s))({\zeta }_t-m(t){)}^T]$, covariance function$t$ 只是一个 index，并不像是 $f(x)$ 中的 $x$ 那样的自变量

1.3 Gaussian Process Regression 高斯过程回归

weight-space view: 关注的是 $w$ $P(y^{\ast }|Data,X^{\ast })={\int }_{w}P(y^{\ast }|w,Data,X^{\ast })P(w)dw$function-space view: 关注的是 $f(x)$ $P(y^{\ast }|Data,X^{\ast })={\int }_{w}P(y^{\ast }|f,Data,X^{\ast })P(f)df$

2 Weight-space view

回顾贝叶斯线性回归（见参考资料 2）

Inference: 求后验 $P(w|Data)=N(w|{\mu }_w,{\Sigma }_w)$ where, $\{\begin{array}{l}{\mu }_w={\sigma }^{-2}{A}^{-1}{X}^{T}Y\\ {\Sigma }_w=A^{-1}\\ A={\sigma }^{-2}{X}^{T}X+{\Sigma }_P^{-1}\end{array}$Prediction: $\text{Given }{x}^{\ast }$ $$\begin{gathered} P(f(x^{\ast })|Data,x^{\ast })=N((x^{\ast }{)}^T{\mu }_w,(x^{\ast }{)}^T{\Sigma }_w{x}^{\ast }) \\ P(y^{\ast }|Data,x^{\ast })=N((x^{\ast }{)}^T{\mu }_w,(x^{\ast }{)}^T{\Sigma }_w{x}^{\ast }+{\sigma }^2) \end{gathered}$$

但如果遇到的模型是非线性的，一般的做法可以是先做一个非线性转换 $z=\varphi (x),x\in {\mathbb{R}}^p,z\in {\mathbb{R}}^q$，然后再使用贝叶斯线性回归（该非线性转换一般把数据由低维转换成高维）

2.1 Process

考虑不存在噪音的情况： $$f(x^{\ast })|X,Y,x^{\ast }\sim N(x^{\ast }({\sigma }^{-2}{A}^{-1}{X}^{T}Y),(x^{\ast }{)}^T{A}^{-1}{x}^{\ast })$$ where, $A={\sigma }^{-2}{X}^{T}X+{\Sigma }_P^{-1}$

$$\begin{gathered} \text{If }\varphi :x\mapsto z,x\in {\mathbb{R}}^p,z=\varphi (x)\in {\mathbb{R}}^q,q>p \\ \text{Define }\Phi =\varphi (X)=[(\varphi (x_1),...,\varphi (x_N){)}^T{]}_{N\times q},X=[(x_1,...,x_N{)}^T{]}_{N\times p},Y=(y_1,...,y_N{)}^T \\ \text{then }f(x)=\varphi (x{)}^{T}w \end{gathered}$$, let ${\varphi }_{\ast }$ denotes $\varphi (x^{\ast })$ $$f(x^{\ast })|X,Y,x^{\ast }\sim N({\varphi }_{\ast }({\sigma }^{-2}{A}^{-1}{\Phi }^{T}Y),{\varphi }_{\ast }^T{A}^{-1}{\varphi }_{\ast })$$ where, $A={\sigma }^{-2}{\Phi }^T\Phi +{\Sigma }_P^{-1}$

如何计算 $A^{-1}$? 使用：Woodbury formula $$(A+UCV{)}^{-1}=A^{-1}-A^{-1}U(C^{-1}+V{A}^{-1}U{)}^{-1}V{A}^{-1}$$

由此可得,

mean: ${\varphi }_{\ast }({\sigma }^{-2}{A}^{-1}{\Phi }^{T}Y)={\varphi }_{\ast }{\Sigma }_p{\Phi }^T(K+{\sigma }^{2}I{)}^{-1}Y$covariance: ${\varphi }_{\ast }^T{A}^{-1}{\varphi }_{\ast }={\varphi }_{\ast }^T{\Sigma }_p{\varphi }_{\ast }-{\varphi }_{\ast }^T{\Sigma }_p{\Phi }^T(K+{\sigma }^{2}I{)}^{-1}\Phi {\Sigma }_p{\varphi }_{\ast }$

where, $K={\Phi }^T{\Sigma }_p\Phi$

2.2 kernel funciton 的意义

观察 $\kappa (x,x^{\prime })$ $\begin{array}{rl}\kappa (x,x^{\prime })& =\varphi (x{)}^T{\Sigma }_p\varphi (x^{\prime })=\varphi (x{)}^T{\Sigma }_p^{\frac{1}{2}}{\Sigma }_p^{\frac{1}{2}}\varphi (x^{\prime })\\ & =({\Sigma }_p^{\frac{1}{2}}\varphi (x){)}^T({\Sigma }_p^{\frac{1}{2}}\varphi (x^{\prime }))\\ & =<\psi (x),\psi (x^{\prime })>\end{array}$，

可见它可以被简化成两数内积的形式。 此时若进一步观察均值和协方差的表达式，可以发现除了 $Y$ 以及 ${\sigma }^{2}I$，其余的所有部分都可以用 kernel function $\kappa (x,x^{\prime })$ 来表达

由此可得，该 kernel 已经包含了原特征空间的所有信息，从而可以成为一种原来方法的简单替代。GP的过程中 kernel function 的重要性就体现于它能将一个复杂问题转化成一个只需要考虑 kernel 的问题

2.3 总结

weight-space view of GPR = Bayesian Linear Regression + Kernel trick (Non-linear Transformation, innner product)

3 From wsv to fsv

先回顾对高斯过程的定义

对于随机变量 $\{{\zeta }_t{\}}_{t\in T},T:\text{continuous time/ space}$

如果对于任意 $\forall n\in {\mathbb{N}}^{+}$, 都存在映射 $t_1,...,t_n\to {\zeta }_1,...,{\zeta }_n$ such that $\mathbf{\zeta }=({\zeta }_1,...,{\zeta }_n{)}^T\sim N(\mathbf{\mu },\mathbf{\Sigma })$, 那么就把 $\{{\zeta }_t{\}}_{t\in T}$ 记为高斯过程, 即 $${\zeta }_t\sim GP(m(t),\kappa (t,s))$$

再回顾 weight-space view（其关注的对象为 $w$）

$\{\begin{array}{l}f(x)=\varphi (x{)}^{T}w\\ y=f(x)+\epsilon ,\epsilon \sim N(0,\sigma )\end{array}$

结合 Section 2.3 中的结论，可以发现 GPR 的 wsv 和 GP 的定义似乎没啥关系。

的确如此，wsv 本质上是贝叶斯线性推断和核方法的结合物，只是能够恰好从另一方面体现高斯过程。而 function-space view 关注的对象则是 $y$，正如 GP 所定义的那样

就具体问题而言，使用 wsv 还是 fsv 没有差别

4 Function-space view

假设样本集 $\{f(x)\}$ 满足 $$\{f(x){\}}_{x\in {\mathbb{R}}^p}\sim GP(m(x),\kappa (x,x^{\prime }))$$

where,

$m(x)=E[f(x)]$$\kappa (x,x^{\prime })=E[(f(x)-m(x))(f(x^{\prime })-m(x^{\prime }){)}^T]$

4.1 Regression

$$\begin{gathered} \text{Data}:\{(x_i,y_i){\}}_{i=1}^N,y=f(x)+\epsilon \\ X=[(x_1,...,x_N{)}^T{]}_{N\times p} \\ Y=[(y_1,...,y_N{)}^T{]}_{N\times 1} \end{gathered}$$

$$\begin{gathered} f(X)\sim N(\mu (X),\kappa (X,X)) \\ Y=f(X)+\epsilon \sim N(\mu (X),\kappa (X,X)+{\sigma }^{2}I) \end{gathered}$$

4.2 Prediiction

$$\begin{gathered} \text{Given}{X}^{\ast }=(x_1^{\ast },...,x_M^{\ast }{)}^T \\ 求Y^{\ast } \\ {Y}^{\ast }=f(X^{\ast })+\epsilon \\ \left(\begin{array}{c}Y\\ f(X^{\ast })\end{array}\right)\sim N\left(\begin{array}{c}\left(\begin{array}{c}\mu (X)\\ \mu (X^{\ast })\end{array}\right),\left(\begin{array}{cc}\kappa (X,X)+{\sigma }^{2}I& \kappa (X,X^{\ast })\\ \kappa (X^{\ast },X)& \kappa (X^{\ast },X^{\ast })\end{array}\right)\end{array}\right) \end{gathered}$$

已知公式，对于 $X\sim N(\mu ,\Sigma )$ $$\begin{gathered} \text{If}X=\left(\begin{array}{c}{X}_a\\ X_b\end{array}\right),\mu =\left(\begin{array}{c}{\mu }_a\\ {\mu }_b\end{array}\right),\Sigma =\left(\begin{array}{cc}{\Sigma }_{aa}& {\Sigma }_{ab}\\ {\Sigma }_{ba}& {\Sigma }_{bb}\end{array}\right) \\ \text{then we have }{X}_b|X_a\sim N({\mu }_{b|a},{\Sigma }_{b|a}) \end{gathered}$$ where,

${\mu }_{b|a}={\Sigma }_{ba}{\Sigma }_{aa}^{-1}(X_a-{\mu }_a)+{\mu }_b$${\Sigma }_{b|a}={\Sigma }_{bb}-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{\Sigma }_{ab}$

根据该公式可得 $$P(f(X^{\ast })|Y,X,X^{\ast })=P(f(X^{\ast })|Y)\sim N({\mu }^{\ast },{\Sigma }^{\ast })$$ where,

${\mu }^{\ast }=\kappa (X^{\ast },X)(\kappa (X,X)+{\sigma }^{2}I{)}^{-1}(Y-\mu (X))+\mu (X^{\ast })$${\Sigma }^{\ast }=\kappa (X^{\ast },X^{\ast })-\kappa (X^{\ast },X)(\kappa (X,X)+{\sigma }^{2}I{)}^{-1}\kappa (X,X^{\ast })$

因此， $$P(y^{\ast }|Y,X,X^{\ast })\sim N({\mu }^{\ast },{\Sigma }^{\ast }+{\sigma }^{2}I)$$

5 Example in code

简单高斯过程回归实现 (斜体部分内容来自参考资料 4)

考虑代码实现一个高斯过程回归，API 接口风格采用 sciki-learn fit-predict 风格。由于高斯过程回归是一种非参数化 (non-parametric) 的模型，每次的 inference 都需要利用所有的训练数据进行计算得到结果，因此并没有一个显式的训练模型参数的过程，所以 fit 方法只需要将训练数据保存下来，实际的 inference 在 predict 方法中进行。

结果如下图，红点是训练数据，蓝线是预测值，浅蓝色区域是 95% 置信区间。真实的函数是一个 cosine 函数，可以看到在训练数据点较为密集的地方，模型预测的不确定性较低，而在训练数据点比较稀疏的区域，模型预测不确定性较高。

python 实现

详见 参考资料 4

matlab 实现

定义类

classdef GPR properties is_fit = 0; train_X = nan; train_y = nan; params = struct('l', 0.5, 'sigma_f', 0.2); end methods function obj = GPR() end function obj = fit(obj, X, y) obj.train_X = X; obj.train_y = y; obj.is_fit = 1; end function [mu, cov] = predict(obj, X) if obj.is_fit == 0 disp('model not fit yet') end Kff = obj.kernel(obj.train_X, obj.train_X); Kyy = obj.kernel(X, X); Kfy = obj.kernel(obj.train_X, X); Kff_inv = inv(Kff + 1e-8 * eye(length(obj.train_X))); mu = Kfy' * Kff_inv * obj.train_y; cov = Kyy - Kfy' * Kff_inv * Kfy; end function result = kernel(obj, x1, x2) distance = x1.^2 + (x2.^2)' - 2 * x1 * x2'; sigma_f = obj.params.('sigma_f'); l = obj.params.('l'); result = sigma_f^2 * exp(-distance/(2*l^2)); end end end

训练 + 预测

train_X = [3, 1, 4, 5, 9]'; train_y = cos(train_X); test_X = 0:0.1:10; test_X = test_X'; gpr = GPR(); gpr = gpr.fit(train_X, train_y); [mu, cov] = gpr.predict(test_X); uncertainty = 1.96 * sqrt(diag(cov)); figure; xfill = [test_X', fliplr(test_X')]; yfill = [(mu - uncertainty)', fliplr((mu + uncertainty)')]; fill(xfill, yfill, 'c', 'FaceAlpha', 0.5, 'EdgeAlpha', 1, 'EdgeColor', 'c'); hold on; plot(test_X, mu);

结果示意图