黑盒优化：从线性回归到高斯过程回归

从线性回归推导高斯过程回归 Rasmussen & Williams (2006)。

1. 线性回归问题

问题描述：给定数据集 $D = \left\{(\mathbf{x}_1, y_1), (\mathbf{x}_2, y_2), ..., (\mathbf{x}_n, y_n)\right\}$，其中 $\mathbf{x}_i = [x_{i1}; x_{i2}; ...; x_{ip}]$ 为 $p$ 维的输入，每维表示一个特征，$y_i \in \mathbf{R}$ 为观测输出。线性回归（linear regression）试图学得一个线性模型，使得输入 $\mathbf{x}_i$ 时，模型的输出能够尽可能的接近观测输出 $y_i$。

通常定义这个线性模型为 $f(\mathbf{x}_i) = \mathbf{w}^T \mathbf{x}_i$，其中 $\mathbf{w}$ 为 $p \times 1$ 的向量。并假设真实观测值与模型输出值相差一个误差项 $\mathbf \epsilon$，即 $y_i = f(\mathbf{x}_i) + \epsilon_i$。进一步假设每一个误差项独立同分布，服从均值为 0，方差为 $\sigma_n^2$ 的高斯分布，即：$\epsilon \sim \mathcal{N}(0, \sigma_n^2)$。

可以通过给 $\mathbf{x}_i$ 添加一个恒为 1 的特征，实现 $f(\mathbf{x}_i) = \mathbf{w}^T \mathbf{x}_i + b$。这里设定 $\mathbf{x}_i$ 的 $p$ 维中已经包含一个恒为 1 的项。

对于包含 $n$ 个数据的数据集 $D$，定义 $\mathbf{X} = [\mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_n]$ 为 $p \times n$ 的矩阵，表示输入矩阵；$\mathbf{y} = [y_1; y_2; ..., y_n;]$ 为 $n \times 1$ 的向量，表示真实输出向量；$\mathbf{\epsilon} = [\epsilon_1; \epsilon_2; ..., \epsilon_n;]$ 为 $n \times 1$ 的向量，表示误差向量。我们可以将模型写成矩阵的形式：

$$f(\mathbf{X}) = \mathbf{X}^T \mathbf{w} \ \ \ \ \mathbf{y} = f(\mathbf{X}) + \mathbf{\epsilon}$$

2. 频率线性回归方法

频率方法完全依靠数据来计算模型的参数 $\mathbf{w}$。最常用的方法是最小二乘法，求一个使得均方误差最小的 $\mathbf{w}^*$：

$$\mathbf{w}^* = \arg\min_{\mathbf{w}} (\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\mathbf{y} - \mathbf{X}^T \mathbf{w})$$

令 $E_\mathbf{w} = (\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\mathbf{y} - \mathbf{X}^T \mathbf{w})$，对 $\mathbf{w}$ 求导得：$\frac{\partial{E_\mathbf{w}}}{\partial{\mathbf{w}}} = 2 \mathbf{X} (\mathbf{X}^T\mathbf{w} - \mathbf{y})$。

令上式为零可以得到 $\mathbf{w}$ 最优解。

$(\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\mathbf{y} - \mathbf{X}^T \mathbf{w}) = \| \mathbf{y} - \mathbf{X}^T \mathbf{w} \|^2$，所以不求导也很显然当 $\mathbf{y} = \mathbf{X}^T \mathbf{w}$ 时，取得最优解。此时

当 $\mathbf{X}^T\mathbf{X}$ 的逆存在的时候，$\mathbf{w}^*$ 可求得唯一解：

$$\begin{aligned} \mathbf{X} \mathbf{X}^T \mathbf{w}^* &= \mathbf{X} \mathbf{y} \\ \mathbf{w}^* &= (\mathbf{X} \mathbf{X}^T)^{-1} \mathbf{X} \mathbf{y} \end{aligned}$$

此时学得的线性回归模型为 $f(\mathbf{x_i}) = \mathbf{x_i}^T (\mathbf{X} \mathbf{X}^T)^{-1} \mathbf{X} \mathbf{y}$。

在问题描述中，假设了真实输出与模型预测相差一个正态分布的误差项，此时最小二乘法求得的最优参数也是其最大似然估计的结果：

当观测数据相互独立时，似然函数为：

$$\begin{aligned} p(\mathbf{y} | \mathbf{X}, \mathbf{w}) & = \prod_{i=1}^{n} p (y_i | \mathbf{x}_i, \mathbf{w}) = \prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi} \sigma_n} \exp \left( - \frac{(y_i - \mathbf{x_i}^T \mathbf{w})^2}{2 \sigma_n^2} \right) \\ &= \frac{1}{(2\pi\sigma_n^2)^{\frac{n}{2}}} \exp \left( - \frac{1}{2 \sigma_n^2} |\mathbf{y} - \mathbf{X}^T \mathbf{w}|^2 \right) \\ &= \frac{1}{(2\pi\sigma_n^2)^{\frac{n}{2}}} \exp \left[ - \frac{1}{2 \sigma_n^2} (\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\mathbf{y} - \mathbf{X}^T \mathbf{w}) \right] \end{aligned}$$

对 $\mathbf{w}$ 求导：$\frac{\partial p(\mathbf{y} | \mathbf{X}, \mathbf{w})}{\partial \mathbf{w}} = p(\mathbf{y} | \mathbf{X}, \mathbf{w}) \times \frac{\partial{E_\mathbf{w}}}{\partial{\mathbf{w}}}$，因为 $p(\mathbf{y} | \mathbf{X}, \mathbf{w}) > 0$，因此令 $\frac{\partial p(\mathbf{y} | \mathbf{X}, \mathbf{w})}{\partial \mathbf{w}} = 0$ 的解等价于最小二乘法中 $\frac{\partial{E_\mathbf{w}}}{\partial{\mathbf{w}}} = 0$ 的解。所以在高斯随机误差的假设下，最小二乘法估计与最大似然估计等价。

小结：

• 频率的线性回归方法中得到的仅仅是模型参数的点估计，预测时使用此点估计值做预测。
• 模型完全由数据决定，即模型的所有信息都编码在可用的训练集中。

3. 贝叶斯线性回归方法

使用贝叶斯方法进行线性回归。首先要给参数指定一个先验分布，其含义为在获取数据之前所掌握的参数信息；然后使用获取的数据依据贝叶斯理论更新掌握的信息，得到参数的后验分布；最后使用此后验分布来做预测。其与频率回归方法相比，有如下两点不同：

• 参数信息不完全由数据决定，还受先验的影响。
• 使用参数的分布进行预测。

具体流程如下：首先指定参数的先验，为了使后验分布能有解析解，这里我们指定似然函数的共轭先验。似然函数是高斯分布，高斯分布与其自身共轭（自共轭），所以指定参数 $\mathbf{w}$ 先验分布为均值为 $\mathbf{0}$，方差为 $\Sigma_{p}$ 的多元高斯分布：$\mathbf{w} \sim \mathcal{N} (\mathbf{0}, \Sigma_{p})$。

参数 $\mathbf{w} \sim \mathcal{N}(\mathbf{0}, \Sigma_p)$ 的先验均值设为 $\mathbf{0}$ 只是为了简单起见，并不一定非得为 $\mathbf{0}$。同时在实际使用中，此设置不会影响实际效果。

获取到数据 $D = (\mathbf{X}, \mathbf{y})$ 之后，数据之间相互独立，且服从正态分布，因此似然函数为：

$$\begin{aligned} p(\mathbf{y} | \mathbf{X}, \mathbf{w}) & = \prod_{i=1}^{n} p (y_i | \mathbf{x}_i, \mathbf{w}) = \prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi} \sigma_n} \exp \left( - \frac{(y_i - x_i^T \mathbf{w})^2}{2 \sigma_n^2} \right) \\ &= \frac{1}{(2\pi\sigma_n^2)^{\frac{n}{2}}} \exp \left( - \frac{1}{2 \sigma_n^2} |\mathbf{y} - \mathbf{X}^T \mathbf{w}|^2 \right) \\ &= \frac{1}{(2\pi\sigma_n^2)^{\frac{n}{2}}} \exp \left[ - \frac{1}{2 \sigma_n^2} (\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\mathbf{y} - \mathbf{X}^T \mathbf{w}) \right] \\ &= \frac{1}{(2\pi\sigma_n^2)^{\frac{n}{2}}} \exp \left[ - \frac{1}{2} (\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\sigma_n^2 \mathit{I})^{-1} (\mathbf{y} - \mathbf{X}^T \mathbf{w}) \right] \\ & \sim \mathcal{N}(\mathbf{X}^T \mathbf{w}, \sigma_n^2 \mathbf{I}) \end{aligned}$$

依据贝叶斯理论：

$$\text{posterior} = \frac{\text{likelihood} \times \text{prior}}{\text{marginal likelihood}}, \quad p(\mathbf{w} | \mathbf{y}, \mathbf{X}) = \frac{p(\mathbf{y} | \mathbf{X}, \mathbf{w}) p(\mathbf{w})}{p(\mathbf{y} | \mathbf{X})}$$

其中 $p(\mathbf{y} | \mathbf{X})$ 是边缘似然函数：

$$p(\mathbf{y} | X) = \int p(\mathbf{y} | X, \mathbf{w}) p(\mathbf{w}) d \mathbf{w}$$

可以推导出参数 $\mathbf{w}$ 的后验分布，此分布包含了先验和似然的所有信息：

$$\begin{aligned} p(\mathbf{w} | \mathbf{y}, \mathbf{X}) = \ & \frac{ p(\mathbf{y} | \mathbf{X}, \mathbf{w}) p(\mathbf{w})}{p(\mathbf{y}|\mathbf{X})} \\ \propto \ & \exp \left( - \frac{1}{2 \sigma_n^2} (\mathbf{y} - \mathbf{X}^T \mathbf{w})^T (\mathbf{y} - \mathbf{X}^T \mathbf{w}) \right) \exp \left( - \frac{1}{2} \mathbf{w}^T \Sigma_p^{-1} \mathbf{w} \right) \\ \propto \ & \exp \left( - \frac{1}{2} (\mathbf{w} - \mathbf{\bar{w}})^T \left(\frac{1}{\sigma_n^2}\mathbf{X}\mathbf{X}^T + \Sigma_{p}^{-1}\right)(\mathbf{w} - \mathbf{\bar{w}}) \right) \\ & \text{其中：} \mathbf{\bar{w}} = \sigma_{n}^{-2}(\sigma_n^{-2}\mathbf{X}\mathbf{X}^T + \Sigma_p^{-1})\mathbf{X}\mathbf{y} \\ \sim \ & \mathcal{N}\left(\mathbf{\bar{w}} = \frac{1}{\sigma_n^2}\mathbf{A}^{-1}\mathbf{X}\mathbf{y}, \ \mathbf{A}^{-1}\right) \\ & \text{其中：} \mathbf{A} = \sigma_{n}^{-2} \mathbf{X} \mathbf{X}^T + \Sigma_p^{-1} \end{aligned}$$

对于一个测试样本 $\mathbf{x}_*$，使用后验预测分布对其输出 $f_{*}$ 进行预测：

$$\begin{aligned} p(f_{*}|\mathbf{x}_{*}, \mathbf{X}, \mathbf{y}) = \ & \int p(f_{*} | \mathbf{x}_{*}, \mathbf{w}) p(\mathbf{w} | \mathbf{X}, \mathbf{y}) d\mathbf{w} \\ = \ & \mathcal{N} \left(\frac{1}{\sigma_n^2} \mathbf{x}_{*}^T \mathbf{A}^{-1} \mathbf{X} \mathbf{y}, \ \mathbf{x}_{*}^T \mathbf{A}^{-1} \mathbf{x}_{*}\right) \end{aligned}$$

后验预测分布仍然是一个高斯分布。可以看出，后验预测分布是对所有参数的预测输出的加权平均（连续情况下是积分的形式）。它不再是使用通过某种策略选择出一个单一的参数进行预测，而是考虑了参数的不确定性（分布）。并且预测的结果也是一个分布，而不是点。 这是与频率方法（点估计）的主要区别。

3.1 核技巧解决数据的非线性

标准贝叶斯线性回归只能表达线性关系，无法捕捉复杂非线性模式。此时可以通过非线性映射将输入空间映射到更高维的特征空间，然后在高维特征空间应用线性模型。只要映射函数和参数 $\mathbf{w}$ 无关，整个模型就仍然还是线性模型（关于参数 $\mathbf{w}$ 是线性的）。另外映射函数如何选？实际上大部分时候不需要显式的做映射，而是使用核技巧绕过。核技巧参考前文「核函数」一节中的讨论。

具体来说，先使用非线性映射函数 $\phi(\mathbf{x})$ 将 $p$ 维的原始输入 $\mathbf{x}$ 映射到一个 $q$ 维的空间中，在新的特征空间中进行线性回归。此时 $\mathbf{w}$ 为 $q \times 1$ 的向量，$\epsilon \sim \mathcal{N}(0, \sigma_n^2)$，回归模型为：

$$f(\mathbf{X}) = \phi(\mathbf{X})^T \mathbf{w}, \quad \mathbf{y} = f(\mathbf{X}) + \mathbf{\epsilon}$$

此时得到后验预测分布为：

$$p(f_{*}|\mathbf{x}_{*}, \mathbf{X}, \mathbf{y}) = \mathcal{N} \left(\frac{1}{\sigma_n^2} \phi(\mathbf{x}_{*})^T \mathbf{A}^{-1} \Phi \mathbf{y}, \quad \phi(\mathbf{x}_{*})^T \mathbf{A}^{-1} \phi(\mathbf{x}_{*})\right)$$

其中 $\Phi = \Phi(\mathbf{X})$，$\mathbf{A} = \sigma_{n}^{-2} \Phi \Phi^T + \Sigma_q^{-1}$。

上面的公式很难实际应用，因为要计算 $\mathbf{A}^{-1}$（$q \times q$ 维），当 $q$ 很大的时候，计算 $\mathbf{A}$ 的逆是一个成本很高的操作，或者当特征空间无限维的时候，根本就无法计算。此时核技巧派上了用场，不仅避免了计算 $\mathbf{A}^{-1}$，连映射函数都无需显式指定了。核技巧是在输入空间直接计算特征空间的内积，因此我们要重写上面的公式，转换成内积的形式。

对于均值：

$$\mathbf{A} \Sigma_q \Phi = \sigma_n^{-2} \Phi\left(\Phi^T \Sigma_q \Phi + \sigma_n^2 I\right)$$

从左侧乘以 $\mathbf{A}^{-1}$，右侧乘以 $\left(\Phi^T \Sigma_q \Phi+\sigma_n^2 I\right)^{-1}$，得到转换前后的均值相等：

$$\Sigma_q \Phi \left(\Phi^T \Sigma_q \Phi+\sigma_n^2 I\right)^{-1} = \sigma_n^{-2} \mathbf{A}^{-1} \Phi$$

对于方差，根据 Woodbury 矩阵恒等式：

$$(\mathbf{Z}+\mathbf{U}\mathbf{W}\mathbf{V}^T)^{-1}=\mathbf{Z}^{-1}-\mathbf{Z}^{-1}\mathbf{U}(\mathbf{W}^{-1}+\mathbf{V}^T\mathbf{Z}^{-1}\mathbf{U})^{-1}\mathbf{V}^T\mathbf{Z}^{-1}$$

令：

$$\mathbf{Z}^{-1} = \Sigma_q, \ \ \ \ \mathbf{W}^{-1} = \sigma_n^2 \mathbf{I}, \ \ \ \ U = V = \Phi$$

可得：

$$(\Sigma_q^{-1} + \sigma_{n}^{-2} \Phi \Phi^T)^{-1} = \Sigma_q - \Sigma_q \Phi (\sigma_n^2 \mathbf{I} + \Phi^T \Sigma_q \Phi)^{-1} \Phi^T \Sigma_q$$

代入可得核技巧转换后后验分布为：

$$\begin{aligned} f_{*}|\mathbf{x}_{*}, \mathbf{X}, \mathbf{y} \sim \mathcal{N} & \left(\phi_{*}^T\Sigma_q \Phi (\mathbf{K} + \sigma^{2}_{n}\mathbf{I})^{-1} \mathbf{y}, \right.\\ & \left. \phi_{*}^{T} \Sigma_q \phi_{*} - \phi_{*}^T \Sigma_q\Phi (\Phi^T \Sigma_q \Phi + \sigma^{2}_{n}\mathbf{I})^{-1} \Phi^T \Sigma_q \phi_{*}\right) \end{aligned}$$

其中 $\phi_* = \phi({\mathbf{x}_*})$，$\mathbf{K} = \Phi^T \Sigma_q \Phi$。这个 $\mathbf{K}$ 刚好是协方差矩阵：

$$\operatorname{Cov}(f(\mathbf{X}), f(\mathbf{X})) = \mathbb{E}\left[f(\mathbf{X}) f\left(\mathbf{X}\right)\right] = \boldsymbol{\Phi}(\mathbf{X})^{\top} \mathbb{E}\left[\mathbf{w} \mathbf{w}^{\top}\right] \boldsymbol{\Phi}\left(\mathbf{X}\right) = \boldsymbol{\Phi^T} \Sigma_{q} \boldsymbol{\Phi}$$

通常情况下测试集包含多个数据，依据上式及协方差的定义，可以将预测写成向量的形式：

$$\begin{aligned} \mathbf{f}_{*} | \mathbf{X}_{*}, \mathbf{X}, \mathbf{y} \sim \mathcal{N} & \left(K\left(\mathbf{X}_{*}, \mathbf{X}\right) \left[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I} \right]^{-1} \mathbf{y}, \right.\\ & \left. K\left(\mathbf{X}_{*}, \mathbf{X}_{*}\right) - K\left(\mathbf{X}_{*}, \mathbf{X}\right) \left[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I} \right]^{-1} K\left(\mathbf{X}, \mathbf{X}_{*}\right)\right) \end{aligned}$$

一些说明：

• 对于每一个测试点 $\mathbf{x}_*$，其后验的均值是 $\mathbf{y}$ 的线性组合；也可以看成是 $k(\mathbf{x}_i, \mathbf{x}_\star)$ 的线性组合：$\bar{f}\left(\mathbf{x}_*\right) = \sum_{i=1}^n \alpha_i k\left(\mathbf{x}_i, \mathbf{x}_*\right)$，其中 $\alpha = \left[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I} \right]^{-1} \mathbf{y}$。
• 对于每一个测试点 $\mathbf{x}_*$，其后验的方差与训练集的输出 $\mathbf{y}$ 无关，仅与输入有关。可以把方差看成是先验方差（$K\left(\mathbf{X}_{*}, \mathbf{X}_{*}\right)$）减去从训练集中获取的信息 $K\left(\mathbf{X}_{*}, \mathbf{X}\right) \left[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I} \right]^{-1} K\left(\mathbf{X}, \mathbf{X}_{*}\right)$，观测数据让方差逐渐变小，也符合逻辑。
• 尽管高斯过程对应无限维特征空间，但每个点的后验分布可表示为仅与训练集相关，因此才能实际应用。
• 协方差矩阵是通过核函数来计算的。核函数的种类有很多，如果我们能够深刻理解数据的特性，可以依此来选择合适的核函数。比如如果明确知道在参数空间上的输出是周期性的，那么可以选择周期性核函数。但在实际使用中的大部分情况下，这些信息是很难知道的，此时没有特别有效的方法来选择核函数。在高斯过程回归中，比较常用的核函数是 RBF 核 和 Matérn 核。

4. 高斯过程回归

4.1 理论

前一节（核技巧）推导的结果实际就是高斯过程回归的结果。本节直接利用高斯过程进行推导。

随机过程是随机变量的集合，高斯过程是随机过程的一种。如果一个随机过程满足任意有限子集服从多元高斯分布，那么此随机过程被称为高斯过程。 一个高斯过程由均值函数（mean function: $m(\mathbf{x})$）和协方差函数（covariance function: $k\left(\mathbf{x}, \mathbf{x}^{\prime}\right)$）唯一确定。用如下符号表示均值函数为 $m(\mathbf{x})$、方差函数为 $k\left(\mathbf{x}, \mathbf{x}^{\prime}\right)$ 的高斯过程 $f(\mathbf{x})$（这里 $f$ 不再表示函数）：

$$f(\mathbf{x}) \sim \mathcal{G} \mathcal{P}\left(m(\mathbf{x}), k\left(\mathbf{x}, \mathbf{x}^{\prime}\right)\right)$$

均值函数和协方差函数定义如下：

$$\begin{aligned} m(\mathbf{x}) &= \mathbb{E} \left[ f(\mathbf{x}) \right] \\ k\left(\mathbf{x}, \mathbf{x}^{\prime}\right) &= \mathbb{E}\left[(f(\mathbf{x})-m(\mathbf{x}))\left(f\left(\mathbf{x}^{\prime}\right)-m\left(\mathbf{x}^{\prime}\right)\right)\right] \end{aligned}$$

对于训练集 $D(\mathbf{X}, \mathbf{y})$，使用高斯过程对其输出进行建模，$\mathbf{y}$ 服从多元高斯分布：

$$\mathbf{y} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{K}(\mathbf{X}, \mathbf{X}) + \sigma_{n}^2 \mathbf{I} \right)$$

对于 $n_*$ 个测试数据 $\mathbf{X}_{*}$，我们想预测这些数据的输出 $\mathbf{f}_*$。根据高斯过程的任意有限子集服从多元高斯分布的性质，$\mathbf{y}$ 与 $\mathbf{f}_*$ 合并后还是多元高斯分布：

$$\left[\begin{array}{c}{\mathbf{y}} \\ {\mathbf{f}_{*}}\end{array}\right] \sim \mathcal{N}\left(\mathbf{0},\left[\begin{array}{cc}{K(\mathbf{X}, \mathbf{X}) + \sigma_{n}^2 \mathbf{I}} & {K\left(\mathbf{X}, \mathbf{X}_{*}\right)} \\ {K\left(\mathbf{X}_{*}, \mathbf{X}\right)} & {K\left(\mathbf{X}_{*}, \mathbf{X}_{*}\right)}\end{array}\right]\right)$$

由于联合分布是正态的，可以很容易能够计算出条件分布。具体参考多元正态分布中的条件分布：

$$\begin{aligned} \mathbf{f}_{*} | \mathbf{X}_{*}, \mathbf{X}, \mathbf{y} \sim \mathcal{N} & \left(K\left(\mathbf{X}_{*}, \mathbf{X}\right) [K(\mathbf{X}, \mathbf{X}) + \sigma_n^2\mathbf{I}]^{-1} \mathbf{y}, \right.\\ & \left. K\left(\mathbf{X}_{*}, \mathbf{X}_{*}\right) - K\left(\mathbf{X}_{*}, \mathbf{X}\right) [K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I}]^{-1} K\left(\mathbf{X}, \mathbf{X}_{*}\right)\right) \end{aligned}$$

这个结果和高斯噪声下的基于核技巧的贝叶斯线性回归推导殊途同归。

对数边缘似然（log marginal likelihood, or evidence） 在调超参的时候会用到（见 §5）。边缘似然是似然与先验乘积的积分：

$$p(\mathbf{y} | \mathbf{X}) = \int p(\mathbf{y} | \mathbf{f}, \mathbf{X}) p(\mathbf{f} | \mathbf{X}) d \mathbf{f}$$

高斯过程回归的先验是高斯过程 $f|\mathbf{X} \sim \mathcal{N}(\mathbf{0}, \mathbf{K})$，因此：

$$\begin{aligned} p(\mathbf{f} | \mathbf{X}) & = \frac{\exp \left(-\frac{1}{2}(\mathbf{f}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1}(\mathbf{f}-\boldsymbol{\mu})\right)}{\sqrt{(2 \pi)^k|\boldsymbol{\Sigma}|}} \\ & = \frac{\exp \left(-\frac{1}{2}(\mathbf{f})^T \boldsymbol{\Sigma}^{-1}(\mathbf{f})\right)}{\sqrt{(2 \pi)^k|\boldsymbol{\Sigma}|}} \end{aligned}$$

取对数可得：

$$\log p(\mathbf{f} \mid \mathbf{X}) = -\frac{1}{2} \mathbf{f}^T \mathbf{K}^{-1} \mathbf{f} - \frac{1}{2} \log |\mathbf{K}| - \frac{n}{2} \log 2 \pi$$

而似然是 $\mathbf{y} | \mathbf{f} \sim \mathcal{N}(\mathbf{f}, \sigma_n^2 \mathbf{I}) \rightarrow \mathbf{y} \sim \mathcal{N}(\mathbf{0}, \mathbf{K} + \sigma_n^2 \mathbf{I})$，因此类比于 $\log p(\mathbf{f} | \mathbf{X})$ 可得：

$$\log p(\mathbf{y} \mid \mathbf{X}) = -\frac{1}{2} \mathbf{y}^{\top}\left(\mathbf{K}+\sigma_n^2 \mathbf{I}\right)^{-1} \mathbf{y} - \frac{1}{2} \log \left|\mathbf{K}+\sigma_n^2 \mathbf{I}\right| - \frac{n}{2} \log 2 \pi$$

4.2 实现

在实现时，直接求逆的时间复杂度高（$O(n^3)$），且容易因数值误差导致结果不准确，因此使用 Cholesky 分解和两次三角矩阵解线性方程组来绕过矩阵求逆。虽然复杂度仍为 $O(n^3)$，但是常数因子更优，约为原来的 1/6；且数值稳定性更高：

$$\begin{aligned} & \mathbf{L} = \text{cholesky}(K(\mathbf{X}, \mathbf{X}) + \sigma_n^2\mathbf{I}) \\ & \text{令} \ \ (\mathbf{L} \mathbf{L}^T)^{-1} \mathbf{y} = \boldsymbol{\alpha} \ \rightarrow \ \boldsymbol{\alpha} = \mathbf{L}^T \backslash (\mathbf{L} \backslash \mathbf{y}) \end{aligned}$$

• $\mathbf{L} \backslash \mathbf{y}$ 表示求解 $\mathbf{L} \mathbf{x} = \mathbf{y}$ 的线性方程组。

可得后验均值：

$$\bar{f}_* = \mathbf{k}(\mathbf{x}_*, \mathbf{X}) \boldsymbol{\alpha}$$

对于方差：

$$\begin{aligned} k\left(\mathbf{x}_{*}, \mathbf{X}\right) [K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 \mathbf{I}]^{-1} k\left(\mathbf{X}, \mathbf{x}_{*}\right) & = k\left(\mathbf{x}_{*}, \mathbf{X}\right) (\mathbf{L} \mathbf{L}^T)^{-1} k\left(\mathbf{X}, \mathbf{x}_{*}\right) \\ & = k\left(\mathbf{X}, \mathbf{x}_{*}\right)^T (\mathbf{L}^T)^{-1} \mathbf{L}^{-1} k\left(\mathbf{X}, \mathbf{x}_{*}\right) \\ & = \left(\mathbf{L}^{-1} k\left(\mathbf{X}, \mathbf{x}_{*}\right)\right)^T \mathbf{L}^{-1} k\left(\mathbf{X}, \mathbf{x}_{*}\right) \end{aligned}$$

令 $\mathbf{v} = \mathbf{L}^{-1} k\left(\mathbf{X}, \mathbf{x}_{*}\right)$，同样通过求解线性方程组求得：

$$\mathbf{v} = \mathbf{L} \backslash K\left(\mathbf{X}, \mathbf{x}_{*}\right)$$

可得后验方差：

$$\mathbb{V}\left[f_*\right] = k\left(\mathbf{x}_*, \mathbf{x}_*\right) - \mathbf{v}^{\top} \mathbf{v}$$

及对数边缘似然：

$$\log p(\mathbf{y} \mid \mathbf{X}) = -\frac{1}{2} \mathbf{y}^{\top} \boldsymbol{\alpha} - \sum_i \log \mathbf{L}_{i i} - \frac{n}{2} \log 2 \pi$$

其中：

• $\frac{1}{2} \log \left|\mathbf{K}+\sigma_n^2 \mathbf{I}\right| = \sum_i \log \mathbf{L}_{i i}$ 利用 Cholesky 分解的下三角性质（$\frac{1}{2}\log|\mathbf{L}\mathbf{L}^T| = \sum_i \log \mathbf{L}_{ii}$）。

5. 模型选择与超参数优化

高斯过程回归的理论很完备，但距离将其变成强大的工具，还有一定的距离。因为整个过程还存在很多不确定性，比如核函数如何选，核函数中的超参数如何确定等等。为使高斯过程回归成为强大的实用工具，必须解决这些问题。我们把这些问题统称为模型选择问题（不同超参数类比为不同的模型），或者叫做模型的训练过程。

5.1 贝叶斯模型选择框架

可以基于贝叶斯理论做模型选择。对模型选择进行分级，抽象成下面三类：

• 参数 $\mathbf{w}$：比如线性回归中的参数 $\mathbf{w}$ 或者神经网络中的权重。
• 超参数 $\boldsymbol{\theta}$：比如核函数中的长度尺度参数，或者神经网络中的权重衰减项（weight decay）。
• 结构 $\mathcal{H}_i$：不同的模型结构，比如线性模型、神经网络等。

对于参数 $\mathbf{w}$，依据贝叶斯理论，其后验分布为：

$$p\left(\mathbf{w} \mid \mathbf{y}, \mathbf{X}, \boldsymbol{\theta}, \mathcal{H}_i\right) = \frac{p\left(\mathbf{y} \mid \mathbf{X}, \mathbf{w}, \mathcal{H}_i\right) p\left(\mathbf{w} \mid \boldsymbol{\theta}, \mathcal{H}_i\right)}{p\left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \mathcal{H}_i\right)}$$

其中分母叫边缘似然（marginal likelihood）：

$$p\left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \mathcal{H}_i\right) = \int p\left(\mathbf{y} \mid \mathbf{X}, \mathbf{w}, \mathcal{H}_i\right) p\left(\mathbf{w} \mid \boldsymbol{\theta}, \mathcal{H}_i\right) d \mathbf{w}$$

对于超参数 $\boldsymbol{\theta}$，其后验分布为：

$$p\left(\boldsymbol{\theta} \mid \mathbf{y}, \mathbf{X}, \mathcal{H}_i\right) = \frac{p\left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \mathcal{H}_i\right) p\left(\boldsymbol{\theta} \mid \mathcal{H}_i\right)}{p\left(\mathbf{y} \mid \mathbf{X}, \mathcal{H}_i\right)}$$

边缘似然为：

$$p\left(\mathbf{y} \mid \mathbf{X}, \mathcal{H}_i\right) = \int p\left(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}, \mathcal{H}_i\right) p\left(\boldsymbol{\theta} \mid \mathcal{H}_i\right) d \boldsymbol{\theta}$$

对于模型结构 $\mathbf{w}$，其后验分布为：

$$p\left(\mathcal{H}_i \mid \mathbf{y}, \mathbf{X}\right) = \frac{p\left(\mathbf{y} \mid \mathbf{X}, \mathcal{H}_i\right) p\left(\mathcal{H}_i\right)}{p(\mathbf{y} \mid \mathbf{X})}$$

因为可选的模型结构通常是离散的集合，因此其边缘似然为：

$$p(\mathbf{y} \mid \mathbf{X}) = \sum_i p\left(\mathbf{y} \mid \mathbf{X}, \mathcal{H}_i\right) p\left(\mathcal{H}_i\right)$$

本节基于贝叶斯理论提供了一个有说服力且统一的模型选择框架。但遗憾的是，对于大多数机器学习模型而言，边缘似然的计算（参数空间上的积分）没有解析解，且没有好的近似方法。不过高斯过程回归是一个例外，其有解析解，下节介绍利用此框架对高斯过程回归做模型选择。

5.2 高斯过程回归中的模型选择

高斯过程回归中的模型选择，主要是指对核函数中的超参数进行调优，比如长度尺度（length-scale）参数。基于贝叶斯理论利用最大化对数边缘似然来做超参调优（公式推导见 §4.1 末尾）：

$$\begin{aligned} \log p(\mathbf{y} \mid \mathbf{X}, \boldsymbol \theta) &= -\frac{1}{2} \mathbf{y}^{\top}\left(\mathbf{K}_{\mathbf{f}}+\sigma_n^2 \mathbf{I}\right)^{-1} \mathbf{y} - \frac{1}{2} \log \left|\mathbf{K}_{\mathbf{f}} +\sigma_n^2 \mathbf{I}\right| - \frac{n}{2} \log 2 \pi \\ &= -\frac{1}{2} \mathbf{y}^{\top}\mathbf{K}_{\mathbf{y}}^{-1} \mathbf{y} - \frac{1}{2} \log \left|\mathbf{K}_{\mathbf{y}}\right| - \frac{n}{2} \log 2 \pi \end{aligned}$$

其中每个部分的含义如下：

• $\frac{1}{2} \mathbf{y}^{\top}\mathbf{K}_{\mathbf{y}}^{-1} \mathbf{y}$：数据拟合项。唯一涉及观测目标的项。
• $\frac{1}{2} \log \left|\mathbf{K}_{\mathbf{y}}\right|$：复杂度惩罚项。仅取决于协方差函数和输入数据 $\mathbf{X}$。
• $\frac{n}{2} \log 2 \pi$：归一化常数。

对对数边缘似然求导可得：

$$\begin{aligned} \frac{\partial}{\partial \theta_j} \log p(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta}) & = \frac{1}{2} \mathbf{y}^{\top} \mathbf{K}^{-1} \frac{\partial \mathbf{K}}{\partial \theta_j} \mathbf{K}^{-1} \mathbf{y} - \frac{1}{2} \operatorname{tr}\left(\mathbf{K}^{-1} \frac{\partial \mathbf{K}}{\partial \theta_j}\right) \\ & = \frac{1}{2} \operatorname{tr}\left(\left(\boldsymbol{\alpha} \boldsymbol{\alpha}^{\top}-\mathbf{K}^{-1}\right) \frac{\partial \mathbf{K}}{\partial \theta_j}\right) \text { where } \boldsymbol{\alpha} = \mathbf{K}^{-1} \mathbf{y} \end{aligned}$$

上式的计算复杂度主要取决于对矩阵 $K$ 的求逆操作，对于 $n \times n$ 的正定对称矩阵求逆的标准方法时间复杂度为 $O(n^3)$。对于每个超参数，偏导数的计算时间复杂度为 $O(n^2)$。整体还算高效。