高斯过程回归浅析
最近因为项目需要学习高斯过程回归的知识，也看了一些知乎大佬的科普文和学术论文，但还是感觉有没理解透彻，直到跟着Hildo Bijl 大神的《Gaussian Process Regression Techniques》手推了一遍 心里才有谱了，这里做个简单的分享。
Bijl文章里面的所有图片都有对应的matlab代码，是非常好的学习资料，强烈推荐去他的GitHub主页去逛逛，可以下载到文章和代码。
先从单变量一个数据点的近似(approximating)分析起，变量符号写成 f 。
1、先验分布
之前可能有一些对该变量的经验（也可以是瞎想，单纯地假设），假设这个变量它一般都出现在零附近，大部分时候在[- \lambda_f ， \lambda_f ]之间波动。因为具体是哪个值不清楚，只有一个大概范围和出现频率的印象，那么用 随机 变量 \underline f 来描述( 用下划线来表示该值是一个随机范围，不带下划线表示该值是确定的 )， \underline f 服从高斯分布：
\underline f \sim(f|m,\lambda_f^2)\tag{1}
其中 m 是随机变量的均值， \lambda_f 是随机变量的标准差。
2、观测（传感器）分布
若现在有一个传感器能观测到该变量的数据。
传感器测量值 \hat f 的结果中包含了噪声， \hat f=f+v ， v 为噪声变量。假设噪声随机变量服从均值为零，方差为 \hat \delta_f^2 的高斯分布 \underline v \sim(\upsilon|0,\hat\delta_f^2) 。则传感器测量的先验关系为： \hat{\underline f} = {\underline f} +\underline v 。
因为传感器测量到的是一个具体的数据，可以不需要下划线来表示它的随机性（是确定的deterministically），则
\underline f=\hat f-\underline v\sim(f|\hat f,\hat \delta_f^2)\tag{2}
（噪声的均值是零，方差为 \hat \delta_f^2 ； \hat f 是确定值，均值是 \hat f ，方差是零）。
3、数据融合
如果能把 先验分布 和 测量分布 结合在一起得到一个新的分布，理所应当的这个新的分布会更加准确。这里应用到一个结论：
定理1 对于同一个随机变量 \underline x ，若有两个概率密度函数 f_\underline x^1(x) 、 f_\underline x^2(x)\ ，且这两个概率密度函数独立，则后验密度函数可表示为：
f_\underline x^{posterior}(x)=\frac{f_\underline x^1(x)f_\underline x^2(x)}{\int_{X}f_\underline x ^1(x^{\prime})f_\underline x^2(x^{\prime})dx^{\prime}}
则数据融合后的 后验分布 服从以下分布：
\underline f\sim\frac{N(m,\lambda_f^2)N(\hat f_1,\hat \sigma_f^2)}{\int_{-\infty}^{\infty}N(m,\lambda_f^2)N(\hat f_1,\hat \sigma_f^2)df}\tag{3}
单变量多个数据点的推广 ， f 此时是一个向量， f=[f_1,f_2...,f_n] ，则先验分布服从多维的形式有：
\underline f\sim(m,K)\tag{4}
m是 n\times1 均值矩阵，K是 n\times n 协方差阵。
【PS】一般把先验分布中的均值和协方差分别用符号 m 和 K 来描述，后验分布中的均值和协方差用 \mu 和 \Sigma 来描述。
对于传感器测量值，服从以下分布：
\underline f\sim N(\hat f,\hat \Sigma_{f})\tag{5}
定理2 和定理1类似，对于随机变量 \underline x ，若它的n次独立测满足分布： \underline x\sim N(m_1,K_1) ... \underline x\sim N(m_n,K_n) ，则它的后验分布服从以下分布： f_{\underline x}^{1:n}(x)=N(\mu_{1:n},\Sigma_{1:n}) ，其中：
\Sigma_{1:n}=(K_1^{-1}+...K_n^{-1})^{-1} , \mu_{1:n}=\Sigma_{1:n}(K_1^{-1}m_1+...K_n^{-1}m_n)^{-1}
根据定理2，则单变量多个数据点服从以下分布：
f\sim N(\mu,\Sigma)\tag{6}
\Sigma=(K_1^{-1}+\hat \Sigma_{f_1}^{-1}+...\hat \Sigma_{f_n}^{-1})^{-1} , \mu=\Sigma(K^{-1}m_1+\hat \Sigma_{f_1}^{-1}+...\hat \Sigma_{f_n}^{-1})^{-1}
4、高斯过程回归
回归分析根据一些已知的数据来预测某些未知数据的值，这是高斯过程的一个主要应用。假设有一组输入 X_m （measurement input set），已知它的观测分布 \underline f_m\sim N(f_m|\hat f_m,\hat \Sigma_m) ，如何根据高斯过程计算另一组输入 X_* (trial input set)的后验分布呢？
和1、2、3节中分析的思路一样，先根据经验写出先验分布，这里用联合分布的形式：
\left[ \begin{matrix} \underline f_m \\ \underline f_* \\ \end{matrix} \right] \sim N( \left[ \begin{matrix} f_m \\ f_* \\ \end{matrix} \right] | \left[ \begin{matrix} m_m \\ m_* \\ \end{matrix} \right] , \left[ \begin{matrix} K_{mm}&K_{m*} \\ K_{*m} &K_{**} \\ \end{matrix} \right] ) \tag{7}
同样的，写出观测（传感器）分布：
\left[ \begin{matrix} \underline f_m \\ \underline f_* \\ \end{matrix} \right] \sim N( \left[ \begin{matrix} f_m \\ f_* \\ \end{matrix} \right] | \left[ \begin{matrix} \hat f_m \\ * \\ \end{matrix} \right] , \left[ \begin{matrix} \hat\Sigma_{fm} &* \\ * & \infty \\ \end{matrix} \right] ) \tag{8}
当前传感器只有输入 X_m 的测量值，而不包含任何 X_* 的观测数据，因此把对应 X_* 的协方差写成无穷大。后续的计算中有求逆步骤， * 处的数据对结果没有影响（受无穷大影响，对应项减为零），可以不关心它到底是什么数据关系。
合并先验分布和观测分布，带入公式（6）中，（推导中主要应用到了矩阵求逆引理），得到高斯回归的一般形式：
\left[ \begin{matrix} \underline f_m \\ \underline f_* \\ \end{matrix} \right] \sim N( \left[ \begin{matrix} f_m \\ f_* \\ \end{matrix} \right] | \left[ \begin{matrix} \mu_m \\ \mu_* \\ \end{matrix} \right] , \left[ \begin{matrix} \Sigma_{mm}&\Sigma_{m*} \\ \Sigma_{*m} &\Sigma_{**} \\ \end{matrix} \right] ) \tag{9}
\left[ \begin{matrix} \Sigma_{mm}&\Sigma_{m*} \\ \Sigma_{*m} &\Sigma_{**} \\ \end{matrix} \right] = \left[ \begin{matrix} K_{mm}-K_{mm}(K_{mm}+\hat \Sigma_{fm})^{-1}K_{mm} & K_{m*}-K_{mm}(K_{mm}+\hat \Sigma_{fm})^{-1}K_{m*} \\ K_{*m}-K_{*m}(K_{mm}+\hat \Sigma_{fm})^{-1}K_{mm} &\ K_{**}-K_{*m}(K_{mm}+\hat \Sigma_{fm})^{-1}K_{m*} \\ \end{matrix} \right]
\left[ \begin{matrix} \mu_{m}\\ \mu_{*} \\ \end{matrix} \right] = \left[ \begin{matrix} m_m+K_{mm}(K_{mm}+\hat \Sigma_{fm})^{-1}(\hat f_m-m_m) \\ m_*+K_{*m}(K_{mm}+\hat \Sigma_{fm})^{-1}(\hat f_m-m_m) \\ \end{matrix} \right]
因为矩阵求逆引理有两种形式，公式（9）还有另外一种结果：
\left[ \begin{matrix} \Sigma_{mm}&\Sigma_{m*} \\ \Sigma_{*m} &\Sigma_{**} \\ \end{matrix} \right] = \left[ \begin{matrix} K_{mm}(K_{mm}+\hat \Sigma_{fm})^{-1}\hat \Sigma_{fm} & \hat \Sigma_{fm}(K_{mm}+\hat \Sigma_{fm})^{-1}K_{m*} \\ K_{*m}(K_{mm}+\hat \Sigma_{fm})^{-1} )\hat\Sigma_{fm} &\ K_{**}-K_{*m}(K_{mm}+\hat \Sigma_{fm})^{-1}K_{m*} \\ \end{matrix} \right]
\left[ \begin{matrix} \mu_{m}\\ \mu_{*} \\ \end{matrix} \right] = \left[ \begin{matrix} \Sigma_{mm}(K_{mm}^{-1}m_m+\hat \Sigma_{fm}^{-1}\hat f_m) \\ m_*+K_{*m}(K_{mm}+\hat \Sigma_{fm})^{-1}(\hat f_m-m_m) \end{matrix} \right]
以上就是高斯过程的一般形式，这里就不把具体的推导过程摆上来了，虽然公式冗长显得很吓人，但只要耐心一点理解起来还是很容易的。再次建议有需要的同学去看这篇论文的附录页，公式推导非常详细。
5、高斯过程拓展
从公式（9）来看，预测 X_* 处的分布需要求解 K_{mm} 的逆，而通常 X_m 的数据量是比较庞大的，这会耗费大量时间。为了解决这个问题，可以考虑 稀疏高斯过程回归 。
稀疏高斯过程回归挑选一组稀疏输入集合 X_u （inducing input set），先用测量数据集 X_m 预测稀疏集 X_u 的后验分布 \underline f_u ，然后利用该分布去预测 X_* 的后验分布，相当于用稀疏集 X_u 替代了原来 X_m 的功能 。 X_u 的维数远小于 X_m 的维数，所以计算时间会大大降低。
虽然稀疏高斯过程回归大大降低了预测 X_* 后验分布的时间，但计算 X_u 的后验分布 \underline f_u 时仍需要计算 K_{mm} 的逆，若此时每个测量值点 \underline f_{m_i} 与稀疏值 \underline f_u 之间概率独立，可以用 FITC方法 进一步提高计算 \underline f_u 过程的速度。
以上的高斯过程回归都是离线计算的，通过适当变换，连续插入测量数据 X_{m_i} 逐步更新均值、协方差参数会得到更快的运算速度。
最后谈一点点个人的看法，在控制和机器人研究方向，我之前接触的都是针对机理的建模方法以及基于模型的控制方法。最近翻了些高斯过程回归的论文，有拿来做状态估计的，也有拿来做黑箱系统辨识的，感觉很新颖，如果能把像高斯过程这类的数据分析方法和传统控制相结合，应该可以做出不错的工程应用案例来。
