一般线性模型

在使用线性模型描述输入输出间关系时，可以用下图表示：

$X$作为输入矩阵，$Y$作为输出，其中$\mathcal w$为权重矩阵。如果用线性方程的形式表述出来：

$$ Y = X^T \mathcal w $$

线性回归就是通过这样一个线性模型的输出与实际值进行比较，并用其结果不断的更新权重矩阵$\mathcal w$以求模型输出最大限度的接近实际值。模型输出我们可以将其称为估计值，通常来说，估计值总是会与实际值有一定的差别，也就是误差，除非实际值与模型输入矩阵$\mathbf{X}$严格呈线性关系，否则这个误差不可避免。我们可以使用一个服从高斯分布的随机变量$\epsilon$表示出估计值与实际值间的误差：

$$ \begin{matrix} \mathbf{y} = X^T \mathcal w + \epsilon & & \epsilon \sim \mathcal N(0, \sigma_n^2 I) \end{matrix}$$

既然有了关于误差（或噪声）的假设，那么接下来就可以引入likelihood的概念了。究竟什么是likelihood？通俗地说，likelihood就是在给定参数的条件下的观察值（也就是上一段所说的实际值）所服从的概率分布或满足的概率密度函数[1]，这个概率密度函数可以将其表示为$p(\mathbf{y}|X, \mathcal w)$。我们可以将估计值$X^T\mathcal w$看做是一个常数$^1$，一个常数加上一个服从均值为0，方差为$\sigma_n^2$的高斯分布的随机变量后服从的依旧是高斯分布，方差不变：

$$ \mathbf{y}|X, \mathcal w \sim \mathcal N(X^T \mathcal w, \sigma_n^2 I) $$

贝叶斯线性回归

在贝叶斯估计理论中，我们先假设权重矩阵服$\mathcal w$服从一个先验分布(prior distribution)，而后再依照训练数据与贝叶斯公式写出其后验分布(posterior distribution)。这里我们可以用$p(\mathcal w)$表示先验分布密度函数，用$p(w|\mathbf{y}, X)$表示后验分布密度函数。先验分布是人为设计的，这里直接假设其为一个均值为0，方差为$\Sigma_p$：

$$ \mathcal w \sim \mathcal N(0, \Sigma_p) $$

依据贝叶斯公式，我们可以写出参数矩阵$\mathcal w$后验分布的表达式[1]：

$$ p(w|\mathbf{y}, X) = \frac{p(\mathbf{y}|X, \mathcal w)p(\mathcal w)}{p(\mathbf{y}|X)} $$

公式中的$p(\mathbf{y}|X)$被称为marginal likelihood，可以发现该分布与参数矩阵$\mathcal w$无关，是与之独立的，所以该分布的对于一个给定观测值的关于参数矩阵$\mathcal w$的后验分布来说相当于是一个常数，并不会影响到后验分布的均值与方差，所以实际上与后验分布相关是likelihood与先验分布的乘积：

$$\begin{align} p(w|\mathbf{y}, X) & = v \exp \left(-\frac{1}{2} \left( (\mathbf{y} - X^T \mathcal w)^T \sigma_n^{-2} I (\mathbf{y} - X^T \mathcal w) \right) \right) \exp\left( -\frac{1}{2} \mathcal w^T \Sigma_p^{-1} \mathcal w \right)\\ & = v\exp\left( -\frac{1}{2} (\mathcal w^T \Sigma_p^{-1} w + (X^T \mathcal w)^T \sigma_n^{-2} I X^T w - (X^T \mathcal w)^T \sigma_n^{-2} I \mathbf{y} - \mathbf{y}^T \sigma_n^{-2} I X^T \mathcal w + \mathbf{y}^T \sigma_n^{-2}I \mathbf{y}) \right) \\ & = v\exp\left( -\frac{1}{2} w^T (\Sigma_p^{-1} + X \sigma_n^{-2} I X^T) w - w^T \mathbf{y} \sigma_n^{-2} I X + C \right) \end{align} $$

如果将上式与高斯分布的一般形式[3]进行对比，可以得到如下三个方程（与求解条件高斯分布的方法类似），假设后验分布的方差为$\Sigma$，均值为$\bar{w}$：
$$ \begin{align} & \Sigma^{-1} = (\Sigma_p^{-1} + X \sigma_n^{-2} I X^T)\\ & \Sigma^{-1} \bar{w} = \mathbf{y} \sigma_n^{-2} I X \\& C = \mathbf{y}^T \sigma_n^{-2} I \mathbf{y} \end{align} $$

由此可以求出后验分布$p(w|\mathbf{y}, X) = \mathcal N (\Sigma \mathbf{y} \sigma_n^{-2} I X, \Sigma)$，其中，$\Sigma = (\Sigma_p^{-1} + X \sigma_n^{-2} I X^T)^{-1}$。对于已得到的参数矩阵$w$，可以将其应用到线性模型的预测中，即给定该线性系统一个新的输入$x_{*}$，得到一个预测分布：

$$\begin{align} & y_{*} = x_{*}^T w \\ & y_{*} \sim \mathcal N(x_{*}^T\bar{w}, x_{*}^T \Sigma x_{*})\end{align}$$

在特征空间中进行预测

有时线性回归在低维数据系统内无法获得良好的预测效果，需要人为构建特征空间或利用某些特殊函数将低维数据映射到高维空间。根据Cover定理，高维空间内数据间的线性特性有更高的概率强于低维空间。在文章【ML系列】“非线性”二元分类中，提到了几个常用的将可以将低维空间转换到高维的映射，这里将这类映射统一用$\phi(X)$来表示。当线性系统的输入为$\phi(X)$时，回归方程$f(X)$与$w$所服从的后验分布可以写为$^2$：

$$\begin{align} & f(X) = \phi(X)^T w \\& w \sim \mathcal N(\Sigma_p \phi(X)(K + \sigma_n^{2}I)^{-1} \mathbf{y},\\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \Sigma_p - \Sigma_p\phi(X)(K + \sigma_n^{2} I)^{-1} \phi(X)^T \Sigma_p) \end{align}$$

这里的$K = \phi(X)^T\Sigma_p \phi(X) $。如果给定新的输入$\phi(x_{*})$，与低维空间中类似，可以得到一个预测分布：

$$ \begin{align} & y_{*} \sim \mathcal N(\phi(x_{*})^T \Sigma_p \phi(X)(K + \sigma_n^{2}I)^{-1} \mathbf{y},\\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \phi(x_{*})^T \Sigma_p\phi(x_{*}) - \phi(x_{*})^T \Sigma_p\phi(X)(K + \sigma_n^{2} I)^{-1} \phi(X)^T \Sigma_p \phi(x_{*})) \end{align}$$

高斯过程回归

在文章【ML系列】条件高斯分布（补充）的最后，提到了条件高斯分布的一个非常重要的应用就是已知协方差与历史数据就可以对新的输入矩阵做预测。现在我们知道，实际上，贝叶斯线性回归完成的是相同的功能，其预测分布实质上也是一个条件概率分布，严格来说应当是随机变量$y_{*}| \mathbf{y}$所服从的分布，这与在条件高斯分布中的描述是一致的。让我们将两种预测分布的方差做一对比，结论就十分显然了。如果贝叶斯线性回归不考虑估计误差方差$\sigma_{n}^2$，那么$y_{*}$与$\mathbf{y}$的联合分布的协方差矩阵就可以写为：
$$\left[\begin{matrix} \phi(X)^T \Sigma_p \phi(X) & \phi(X)^T \Sigma_p \phi(x_{*})\\ \phi(x_{*})^T \Sigma_p \phi(X) & \phi(x_{*})^T \Sigma_p \phi(x_{*}) \end{matrix}\right]$$

使用协方差的定义也一样可以得到相同的结果，比如计算预测值$y_{*}$与观测值$\mathbf{y}$的协方差（因为不考虑误差，默认$y_{*} = f(x_{*})$，$\mathbf{y} = f(X)$，因为权重矩阵$w$是服从高斯分布的，所以$f(X)$同服从高斯分布）：

$$ \mathbb{E} (f(x_{*})f(X)^T) = \phi(x_{*})^TE(ww^T)\phi(X) = \phi(x_{*})^T\Sigma_p\phi(X) $$

从这里我们可以发现，由贝叶斯估计得出的结论可以与高斯过程完美融洽。高斯过程回归也就是从高斯过程的角度来理解贝叶斯估计。这样一来，我们可以将高斯联合分布的协方差矩阵视为特征空间映射与先验分布协方差矩阵的内积，但是从定义上来说高斯过程回归是不考虑先验分布的，那么如何来求解预测分布所需要的协方差矩阵呢？这时就可以引入核(kernel)的概念了。这里可以定义一个核函数$k(x, x)$：

$$ k(x, x) = \phi(x)^T\Sigma_p\phi(x) $$

因为先验分布的协方差矩阵$\Sigma_p$是一个实对称矩阵，所以一定可以写为$\Sigma_p = (\Sigma_p^{\frac{1}{2}})^2$。这样的话，就可以将这个核函数拆为两个部分，第一个部分为$\phi(x)^T\Sigma_p^{\frac{1}{2}}$，第二个部分为$\Sigma_p^{\frac{1}{2}}\phi(x)$，这两部分实际上互为转置矩阵，所以可以定义其中一个部分为$\psi(x)$，那么核函数就可以写为：

$$ k(x, x) = \psi(x)^T \psi(x) $$

在经过变换后得到的新的映射$\psi$依旧可以看作是可将低维数据转变到高维的映射，因为从定义上权重矩阵与$ w $与特征空间映射$\phi$的维度是相同的，即使与先验分布协方差进行内积运算，维度也不会发生改变。

在实际应用高斯过程回归时，我们并不会先寻找特征空间映射，然后再计算核函数，而是利用已有的某些特殊的内积映射直接计算核函数，理论上来说，这些特殊的内积映射都可以分解成特征空间映射$\psi$与$\psi^T$的内积。在文章【ML系列】“非线性”二元分类中我提到过两个特殊的内积映射，一是n次多项式内积映射，二是高斯内积映射，特别地，对于高斯内积映射来说，其特征空间映射$\psi$可以将低维空间数据映射到无限维空间中（证明过程日后会补充）。

现在让我们总结一下高斯过程回归的应用步骤：

1. 选取合适的核函数，并计算核矩阵

2. 将核矩阵作为联合高斯分布的协方差矩阵，与历史数据联合计算条件概率分布

3. 由条件概率分布做线性回归预测

下图给出了一个高斯过程回归的应用实例，阴影部分为置信度95%的估计区间，蓝色点为历史数据，红线为预测线

[1] C. E. Rasmussen & C. K. I. Williams, Gaussian Processes for Machine Learning, the MIT Press, 2006, ISBN 026218253X