高斯过程回归（Gaussian Process Regression）

在概率论和统计学中，高斯过程是指观测发生在连续域（例如：时域、空间域）中的一种特殊的概率模型

1 基本概念

在高斯过程，连续的输入空间的任何点与正态分布的随机变量相关，而且任何随机变量的有限集合满足多重正态分布，例如变量间的任意线性组合是正态分布，高斯过程分布是所有随机变量在连续域中的联合分布
在机器学习理论中，针对于推广训练数据（generalize train data）的算法，如果学习方法在对系统发出请求之前进行，称为急切学习（eager learning），如果学习方法滞后于对系统的请求，称为懒惰学习（lazing learning），K近邻算法就是典型的懒惰学习方法
懒惰学习方法一般模型可移植性强，适用于不同的问题域，但要求较大空间来实时保存训练数据集，因为缺少对数据集合的抽象而使得预测过程计算量大。懒惰分类器适用于数据集的特征较少的情况
因此，从机器学习算法的观点看，高斯过程采用懒惰学习的方式来度量样本间的相似度（核函数），从而基于训练数据预测位置观测输入的结果。预测结果往往是边缘分布（多个概率分布函数的累加）
函数空间视角（function-space view）：高斯过程可以看作是基于多个函数定义分布，并在函数空间做推断分析

2 问题描述

对于给定数据集DN={xn,yn}n=1ND_N=\{\mathbf{x}_n,y_n\}_{n=1}^NDN={xn,yn}n=1N，其中xn∈ℜm\mathbf{x}_n\in \Re^mxn∈ℜm，令X=[x1,...,xN]TX=[\mathbf{x}_1,...,\mathbf{x}_N]^TX=[x1,...,xN]T作为观测输入数据矩阵，即输入空间；Y=[y1,…,yN]TY=[y_1,\dots,y_N]^TY=[y1,…,yN]T作为观测输出向量，即输出空间。考虑非线性映射ϕ(x):x∈ℜm→F\phi(\mathbf{x}):\mathbf{x}\in \Re^m \rightarrow Fϕ(x):x∈ℜm→F，这个映射可能是未知的也可能是无限维的。核函数k(xi,xj)k(\mathbf{x}_i,\mathbf{x}_j)k(xi,xj)可以满足这个特性，核函数用内积形式来构建特征空间F为：k(xi,xj)=<ϕ(xi),ϕ(xj)>k(\mathbf{x}_i,\mathbf{x}_j)=<\phi(\mathbf{x}_i),\phi(\mathbf{x}_j)>k(xi,xj)=<ϕ(xi),ϕ(xj)>典型的核函数包括径向基函数k(xi,xj)=exp⁡(−∥xi−xj∥2/2ρ2)k(\mathbf{x}_i,\mathbf{x}_j)=\exp(-\Vert\mathbf{x}_i-\mathbf{x}_j\Vert^2/2\rho^2)k(xi,xj)=exp(−∥xi−xj∥2/2ρ2)，其中ρ>0\rho>0ρ>0是宽度参数
定义高斯分布N(μ,Σ)N(\mu,\Sigma)N(μ,Σ)，其中μ\muμ是均值，Σ\SigmaΣ是协方差
在高斯回归过程模型（Gaussian process regression，GPR）中，每个样本yny_nyn可以表示为：y=f(x)+εy=f(\mathbf{x})+\varepsilony=f(x)+ε其中fff是零均值高斯过程f∼N(0,KXX)\mathbf{f}\sim N(\mathbf{0},K_{XX})f∼N(0,KXX)，KXX={k(xi,xj)}∈ℜN×NK_{XX}=\{k(\mathbf{x}_i,\mathbf{x}_j)\}\in \Re^{N \times N}KXX={k(xi,xj)}∈ℜN×N是核函数定义的特定协方差矩阵，ε∼N(0,σ2)\varepsilon \sim N(0,\sigma^2)ε∼N(0,σ2)
定义KXx=[k(x1,x),…,k(xN,x)]T∈ℜN\mathbf{K}_{X\mathbf{x}}=[k(\mathbf{x}_1,\mathbf{x}),\dots,k(\mathbf{x}_N,\mathbf{x})]^T\in \Re^NKXx=[k(x1,x),…,k(xN,x)]T∈ℜN
典型的高斯过程回归方法是为了对于任意测试数据x∗∈X\mathbf{x}^*\in Xx∗∈X估计预测分布p(y∣x∗)p(y|\mathbf{x}^*)p(y∣x∗)，基于高斯概率假设，可以得到条件分布：p^(y∣x∗,X,Y)∼N(f(x∗),g(x∗))\hat{p}(y|\mathbf{x}^*,X,Y)\sim N(f(\mathbf{x}^*),g(\mathbf{x}^*))p^(y∣x∗,X,Y)∼N(f(x∗),g(x∗))这里，f(x∗)=KXx∗T(KXX+σ2I)−1Yf(\mathbf{x}^*)=\mathbf{K}_{X\mathbf{x}^*}^T(K_{XX}+\sigma^2\mathbf{I})^{-1} Yf(x∗)=KXx∗T(KXX+σ2I)−1Yg(x∗)=σ2+k(x∗,x∗)−kXx∗T(KXX+σ2I)−1kXx∗g(\mathbf{x}^*)=\sigma^2+k(\mathbf{x}^*,\mathbf{x}^*)-\mathbf{k}_{X\mathbf{x}^*}^T(K_{XX}+\sigma^2\mathbf{I})^{-1}\mathbf{k}_{X\mathbf{x}^*}g(x∗)=σ2+k(x∗,x∗)−kXx∗T(KXX+σ2I)−1kXx∗其中I\mathbf{I}I是相应规模的单位矩阵
特别地，令a=[a1,…,aN]T=(KXX+σ2I)−1Y\mathbf{a}=[a_1,\dots,a_N]^T=(K_{XX}+\sigma^2\mathbf{I})^{-1} Ya=[a1,…,aN]T=(KXX+σ2I)−1Y。因此f(x∗)=aTkXx∗=∑i=1Naik(xi,x∗)f(\mathbf{x}^*)=\mathbf{a}^T\mathbf{k}_{X\mathbf{x}^*}=\sum_{i=1}^N a_ik(\mathbf{x}_i,\mathbf{x}^*)f(x∗)=aTkXx∗=∑i=1Naik(xi,x∗)
由此可见，高斯过程可以用一系列的基函数组合来表示

3 高斯过程模型估计

在高斯模型估计中，噪声方差往往视为一个参数，并于核函数在在一起，即σ\sigmaσ
常用的模型参数估计方法是边际概率p(Y∣X)p(Y|X)p(Y∣X)，可以通过条件概率和先验的乘积积分来表示，即为p(Y∣X)=∫p(Y∣f,X)p(f∣X)dfp(Y|X)=\int p(Y|\mathbf{f},X)p(\mathbf{f}|X)d\mathbf{f}p(Y∣X)=∫p(Y∣f,X)p(f∣X)df边际概率取对数有：JML=log⁡p(Y∣X)=−12YT(KXX+σ2I)−1Y−12log⁡det⁡(KXX+σ2I)−N2log⁡(2π)J^{ML}=\log p(Y|X)=-\frac{1}{2}Y^T(\mathbf{K}_{XX}+\sigma^2\mathbf{I})^{-1}Y-\frac{1}{2}\log \det(\mathbf{K}_{XX}+\sigma^2\mathbf{I})-\frac{N}{2}\log(2\pi)JML=logp(Y∣X)=−21YT(KXX+σ2I)−1Y−21logdet(KXX+σ2I)−2Nlog(2π)
或者，考虑输出概率密度函数和高斯过程回归估计之间的Kullbak-Leibler（K-L）散度作为代价函数，即有KL=∫p(y)log⁡p(y)p^(y∣X,Y)dy=∫p(y)log⁡p(y)dy−∫log⁡p^(y∣X,Y)p(y)dyKL=\int p(y)\log\frac{p(y)}{\hat{p}(y|X,Y)}dy=\int p(y)\log p(y)dy-\int\log\hat{p}(y|X,Y)p(y)dyKL=∫p(y)logp^(y∣X,Y)p(y)dy=∫p(y)logp(y)dy−∫logp^(y∣X,Y)p(y)dy其中上式第二项R=∫log⁡p^(y∣X,Y)p(y)dy≈E(log⁡p^(y∣X,Y))求最大值R=\int\log\hat{p}(y|X,Y)p(y)dy\approx E(\log\hat{p}(y|X,Y))求最大值R=∫logp^(y∣X,Y)p(y)dy≈E(logp^(y∣X,Y))求最大值.
由于p^(y∣X,Y)=∫p^(y∣x,X,Y)p(x)d(x)=E(p^(y∣x,X,Y))≈1N∑j=1N12πg(xj)exp⁡(−(y−f(xj))22g(xj))\hat{p}(y|X,Y)=\int \hat{p}(y|\mathbf{x},X,Y)p(\mathbf{x})d(\mathbf{x})=E(\hat{p}(y|\mathbf{x},X,Y))\approx \frac{1}{N}\sum_{j=1}^N \frac{1}{\sqrt{2\pi g(\mathbf{x}_j)}}\exp \left(-\frac{(y-f(\mathbf{x}_j))^2}{2g(\mathbf{x}_j)}\right)p^(y∣X,Y)=∫p^(y∣x,X,Y)p(x)d(x)=E(p^(y∣x,X,Y))≈N1j=1∑N2πg(xj)1exp(−2g(xj)(y−f(xj))2)因此，将上式带入RRR，并引入p(x)p(\mathbf{x})p(x)的真实密度，可以得到R≈JKL=1N∑i=1Nlog⁡(1N∑j=1N12πg(xj)exp⁡(−(y−f(xj))22g(xj)))=1N∑i=1Nlog⁡(1N∑j=1Npi,j)R\approx J^{KL}=\frac{1}{N}\sum_{i=1}^N\log\left(\frac{1}{N}\sum_{j=1}^N \frac{1}{\sqrt{2\pi g(\mathbf{x}_j)}}\exp \left(-\frac{(y-f(\mathbf{x}_j))^2}{2g(\mathbf{x}_j)}\right)\right)=\frac{1}{N}\sum_{i=1}^N\log\left(\frac{1}{N}\sum_{j=1}^Np_{i,j}\right)R≈JKL=N1i=1∑Nlog(N1j=1∑N2πg(xj)1exp(−2g(xj)(y−f(xj))2))=N1i=1∑Nlog(N1j=1∑Npi,j)其中pi,j=12πgjexp⁡(−ei,j22gj)p_{i,j}=\frac{1}{\sqrt{2\pi g_j}}\exp\left(-\frac{e_{i,j}^2}{2g_j}\right)pi,j=2πgj1exp(−2gjei,j2)，这里ei,j=yi−fje_{i,j}=y_i-f_jei,j=yi−fj，fif_ifi和gjg_jgj分别定义f(xi)f(\mathbf{x}_i)f(xi)和g(xj)g(\mathbf{x}_j)g(xj)
由于真实密度p(y)p(y)p(y)和p(x)p(\mathbf{x})p(x)并不知道，但不同级别的概率函数可以调整至结晶误差的N−1/2N^{-1/2}N−1/2阶，所以可以用估计器插值法近似KL散度
高斯过程参数估计不管是通过JKLJ^{KL}JKL还是JMLJ^{ML}JML，都是非线性优化任务，需要对(KXX+σ2I)−1(K_{XX}+\sigma^2\mathbf{I})^{-1}(KXX+σ2I)−1进行迭代，每步计算复杂度是O(N3)O(N^3)O(N3)

4 带噪声的离散观测值估计实例

根据上面内容，已知有先验分布y∼N(0,KXX+σ2I)y\sim N(0,K_{XX}+\sigma^2\mathbf{I})y∼N(0,KXX+σ2I)，对于任意测试数据y∗y^*y∗有：
[yy∗]∼[KXX+σ2IKXxKXxTKxx]\left[ \begin{array}{l} y\\ y^* \end{array} \right] \sim \left[ \begin{array}{l} \mathbf{K}_{XX}+\sigma^2\mathbf{I}&\mathbf{K}_{X\mathbf{x}}\\ \mathbf{K}_{X\mathbf{x}}^T&\mathbf{K}_{\mathbf{x}\mathbf{x}} \end{array} \right][yy∗]∼[KXX+σ2IKXxTKXxKxx]
其中，Kxx\mathbf{K}_{\mathbf{x}\mathbf{x}}Kxx恒等于1

附 matlab写的高斯滤波器

real1 = TestingData_output;
out1 = testingoutput';
N = length(out);
in = zeros(N,1);
in(:,1) = 1:N;
out = (out1 - mean(out1))/std(out1);
real = (real1 - mean(real1))/std(real1);
sigma_square = 0.5;
rho_opt = 100;N = length(in);
K =zeros(N,N);
for i = 1:Nfor j =1:NK(i,j)=exp(-((in(i,1)-in(j,1))^2)/(2*rho_opt^2));end
end
Num = N;
cnt = 0;
K_new = zeros(N,1);
ret = zeros(Num,1);
temp = inv(K+sigma_square*eye(N,N))*out;
for i = 1:Numx =in(i,1);   for k = 1:NK_new(k,1)=exp(-((x-in(k,1))^2)/(2*rho_opt^2));endcnt =cnt + 1;ret(cnt,1) = K_new' * temp;
end
ret =  ret * std(out1) + mean(out1);
plot(ret);hold on
plot(real1); hold on
plot(out1);
%%%%
legend('ret','TestingData_output','testingoutput');