【摘 要】高斯过程作为一种用于预测的非参数模型，可以用于回归任务，也可以用于分类任务。高斯过程面临的最大问题在于，当面临大数据时，其计算复杂度为 O(n3)\mathcal{O}(n^3)，内存复杂度为 O(n2)\mathcal{O}(n^2)，这使其在新形势下的使用非常困难。本文对2006年之前的主要大数据的高斯过程处理方法进行了综述，可以作为了解该方向工作的基础。《机器学习中的高斯过程》一书是高斯过程研究领域的扛鼎之作，本文主要节选自该书的第八章。

【原 文】 Rasmussen, C.E. and Williams, C.K. (2006), Chapter 8 of Gaussian processes for machine learning. Cambridge, Mass: MIT press Cambridge, MA (3).

正如我们在前面的章节中看到的，高斯过程预测的一个重要问题是它的计算规模通常为 \mamthcal{O}(n^3)。对于大型问题（例如 n>10,000n > 10,000），存储 Gram 矩阵和求解相关的线性系统在现代工作站上都是很难的。

已经提出了广泛的建议来处理这个问题。下面我们将它们分为五个部分：

  • 在 第 8.1 节 中，我们考虑 Gram 矩阵的降秩近似；
  • 在 第 8.2 节 中，我们描述了贪婪近似的一般策略；
  • 在 第 8.3 节 中，我们讨论了针对固定超参数近似高斯过程回归问题的各种方法；
  • 在 第 8.4 节 中，我们描述了用于近似固定超参数的高斯过程分类问题的各种方法；
  • 在 第 8.5 节 中，我们描述了近似边缘似然及其导数的方法。

许多这些方法使用训练示例中大小为 m<nm<n 的子集。

8.1 Gram矩阵的降阶近似

在高斯过程回归问题中，我们需要求逆矩阵 K+σn2IK + \sigma^2_nI（或者至少求解线性系统 (K+σn2I)v=yforv(K + \sigma^2_nI)v = y for v）。如果矩阵 KK 的秩为 qq（这样它可以表示为 K=QQ⊤K = QQ^{\top}，其中 QQ 是 n×qn × q 矩阵），则可以使用矩阵求逆引理 式 (A.9) 加速矩阵求逆。作为 (QQ⊤+σn2In)−1=σn−2In−σn−2Q(σn2Iq+Q⊤Q)−1Q⊤(QQ^{\top} + \sigma^2_n I_n)^{-1} = \sigma^{-2}_n I_n − \sigma^{-2}_n Q(\sigma^2_n I_q + Q^{\top} Q)^{-1} Q^{\top}。请注意，n×nn × n 矩阵的求逆现在已转换为 q×qq × q 矩阵的求逆。

如果核是从具有 NN 个特征的显式特征扩展派生的，那么 Gram 矩阵将具有秩 min⁡(n,N)\min(n, N )，因此如果 n>Nn > N ，则利用此结构将是有益的。即使核是非退化的，它也可能具有快速衰减的特征谱（参见例如第 4.3.1 节），因此降秩近似将是准确的。

如果 KK 的秩不 <n< n，我们仍然可以考虑对 KK 进行降秩近似。K 相对于 Frobenius 范数的最优降秩近似 （参见等式 (A.16)）是 UqΛqUq⊤U_q \Lambda_q U^{\top}_q ，其中 Λq\Lambda_q 是 KK 的前导 qq 个特征值的对角矩阵，UqU_q 是相应的正交特征向量的矩阵 [Golub 和 Van Loan，1989，定理8.1.9]。不幸的是，这在实践中的意义有限，因为计算特征分解是一个 \mamthcal{O}(n^3) 操作。然而，它确实表明，如果我们可以更便宜地获得近似特征分解，那么这可能会产生对 KK 的有用的降秩近似。

我们现在考虑从 nn 个数据点中选择一个子集 II；集合 II 的大小为 m<nm<n。剩余的 n−mn − m 个数据点形成集合 RR。（作为助记符，II 代表包含的数据点，RR 代表剩余的点。）我们有时将包含的集合称为活动集。在不失一般性的情况下，我们假设数据点是有序的，因此集合 II 排在第一位。因此 KK 可以划分为

K=[KmmKm(n−m)K(n−m)mK(n−m)(n−m)](8.1)K= \begin{bmatrix} K_{mm} & K_{m(n-m)} \\ K_{(n-m)m} & K_{(n-m)(n-m)} \end{bmatrix} \tag{8.1}

顶部的 m×nm×n 块也称为 KmnK_{mn}，其转置称为 K_{nm}。

在第 4.3.2 节中，我们了解了如何使用 Nystrom 方法来近似核的特征函数。我们现在可以应用相同的想法来近似 KK 的特征值/向量。我们计算 KmmK_{mm} 的特征向量和特征值，并将它们表示为 {λi(m)}i=1m\{\lambda^{(m)}_i \}^{m}_{i=1} 和 {ui(m)}i=1m\{\mathbf{u}^{(m)}_i \}^{m}_{i=1}。使用等式将这些扩展到所有 nn 个点。式 (4.44) 给出

\tide{λ}^{(n)}_i \triangleq \frac{n}{m} \lambda^{(m)}_i, \qquad i = 1, \ldots , m tag{8.2}

u~i(n)≜sqrtmn1λi(m)Knmui(m),i=1,…,m(8.3)\tilde{\mathbf{u}}^{(n)}_i \triangleq sqrt{\frac{m}{n}} \frac{1}{\lambda^{(m)}_i} K_{nm} \mathbf{u}^{(m)}_i, \qquad i = 1, \ldots , m \tag{8.3}

其中已选择 u~i(n)\tilde{\mathbf{u}}^{(n)}_i 的缩放比例，使得 ∣u~i(n)∣≃1|\tilde{\mathbf{u}}^{(n)}_i | \simeq 1。一般来说，我们可以选择在 KK 的近似值中包含多少个近似特征值/向量；选择第一个 pp 我们得到 K~=∑i=1pλ~i(n)u~i(n)(u~i(n))⊤\tilde{K} = \sum^{p}_{i=1} \tilde{\lambda}^{(n)}_i \tilde{\mathbf{u}}^{(n)}_i (\tilde{\mathbf{u}}^{(n)}_i )^{\top}。下面我们将 p=mp = m 得到

K~=KnmKmm−1Kmn(8.4)\tilde{K} = K_{nm}K^{-1}_{mm} K_{mn} \tag{8.4}

使用方程 8.2 和 8.3，我们称之为 K 的 Nystrom 近似。计算 \tilde{K} 需要时间 \mamthcal{O}(m2n)，因为 K_{mm} 的特征分解是 \mamthcal{O}(m3)，每个\tilde{\mathbf{u}}^{(n)}_i 的计算是 \mamthcal{O}(锰）。福克斯等人。 [2001] 已经应用 Nystrom 方法来近似计算机视觉问题中的前几个特征向量，其中所讨论的矩阵的大小大于 106×106。

上面已应用 Nystrom 近似来近似 K 的元素。但是，使用第 i 个特征函数的近似 ̃ φi(x) = (√m/\lambda^{(m)}_i )km(x)^{\top}\mathbf{u}^{(m)}_i , 其中 km(x) = (k(x, x1), . . , k(x, xm))^{\top}（使用当前符号重述等式 (4.44)）和 λi ’ \lambda^{(m)}_i / m 很容易看出，通常我们得到核的近似值 k(x, x′) = ∑N i=1 λiφi(x)φi(x′) 作为i=1 \lambda^{(m)}_i m ̃ φi(x) ̃ φi(x′) (8.5) = m ∑ i=1 \lambda^{(m)}_i m m (\lambda^{(m)}_i )2 km(x)^{\top}\mathbf{u}^{(m)}_i (\mathbf{u}^{(m)}i )^{\top}km(x′) (8.6) = km(x)^{\top}K^{-1}{mm}km(x′)。 (8.7)

显然当量。 (8.4) 是通过评估等式获得的。 (8.7) 对于训练集中的所有数据点对。

通过乘以 eq。 (8.4) 利用 K_{mn} = [K_{mm}K_{m(n-m)}] 容易证明 K_{mm} = \tilde{K}{mm}, K{m(n-m)} = \tilde{K}{m(n-m)}, K{(n-m)m} = \tilde{K} (n−m)m，但是 \tilde{K}{(n-m)(n-m)} = K{(n-m)m}K^{-1}{mm}K{m(n-m)}。区别

K_{(n-m)(n-m)} − \tilde{K}(n−m)(nm) 实际上是 K_{mm} 的 Schur 补集 [Golub and Van Loan, 1989, p. 103]。容易看出 K_{(n-m)(n-m)} − \tilde{K}_{(n-m)(n-m)} 是半正定的；如果向量 f 被划分为 f > = (f > m, f> n−m) 且 f 服从均值为零且协方差为 K 的高斯分布，则 fn−m|fm 将 Schur 补码作为其协方差矩阵，参见等式. （A.6）。

Nystrom 近似是由 Williams 和 Seeger [2001] 以上述方式导出的，用于核机器的应用。 Smola 和 Sch ̈olkopf [2000]（以及 Sch ̈olkopf 和 Smola [2002，第 10.2 节]）提出了另一种产生相同近似值的观点。这里的出发点是我们希望将以点 xi 为中心的核近似为来自活动集的核的线性组合，使得 k(xi, x) ’ ∑ j∈I cijk(xj, x) , ^ k( xi, x) (8.8)

对于一些要确定的系数 {cij} 以优化近似。一个合理的最小化标准是 E© = n ∑ i=1 ‖k(xi, x) − ^ k(xi, x)‖2H (8.9) = tr K − 2 tr(CK_{mn}) + tr(CK_{mm}C^{\top} ), (8.10)

其中系数排列成 n × m 矩阵 C。最小化 E© w.r.t. C 给出 Copt = K_{nm}K^{-1}{mm}；因此我们得到与等式一致的近似值 ^ K= K{nm} K −1 mmK_{mn}。 (8.4)。此外，可以证明 E(Copt) = tr(K − ^ K )。

Smola 和 Sch ̈olkopf [2000] 建议使用贪婪算法来选择要包含在活动集中的点，以最小化误差标准。由于包含一个新数据点需要 \mamthcal{O}(mn) 次操作来评估 E 的变化（参见练习 8.7.2），因此在每次迭代中考虑包含集合 R 的所有成员是不可行的；相反，Smola 和 Sch ̈olkopf [2000] 建议在每次迭代中从集合 R 的随机选择子集中找到最佳点。

Drineas 和 Mahoney [2005] 最近的工作分析了一个与 Nystrom 近似类似的算法，不同之处在于它们使用带替换的偏置采样（选择 K 的第 i 列的概率为 ∝ k2 ii）和内部 m × m 矩阵的伪逆。对于这个算法，他们能够提供近似质量的概率界限。 Frieze 等人的早期工作。 [1998] 使用其行和列的加权随机子采样以及概率误差界限开发了矩形矩阵的奇异值分解 (SVD) 的近似值。然而，这与 Nystrom 近似有很大不同；参见 Drineas 和 Mahoney [2005，sec. 5.2]了解详情。

Fine 和 Scheinberg [2002] 建议使用不完全 Cholesky 分解对 K 进行替代低秩近似（参见 Golub 和 Van Loan [1989，第 10.3.2 节]）。这里的想法是，当计算低于某个阈值的 K 个枢轴的 Cholesky 分解时，将被跳过。2 如果大于阈值的枢轴数为 k，则不完全的 Cholesky 分解需要时间 \mamthcal{O}(nk2)。

8.2 贪心近似

下面描述的许多方法使用从大小为 n > m 的训练集中选择的大小为 m 的活动训练点集。我们假设由于组合学的原因不可能搜索大小为 m 的最优子集。活动集中的点可以随机选择，但一般来说，如果这些点是贪婪地选择 w.r.t.，我们可能会期望更好的性能。一些标准。在统计文献中，贪婪方法也称为前向选择策略。

算法 8.1 中给出了贪婪近似的一般方法。该算法从活动集 I 为空开始，集 R 包含所有训练示例的索引。在每次迭代中，从 R 中选择一个索引并添加到 I。这是通过评估某些标准 Δ 并选择优化该标准的数据点来实现的。对于某些算法，在 R 中的所有点上评估 Δ 可能过于昂贵，因此可以选择一些工作集 J ⊂ R，通常是从 R 中随机选择。贪心选择方法已用于回归子集 (SR)，数据点子集 (SD) 和下面描述的预计过程 (PP) 方法。

8.3 固定超参数的 GPR 近似

我们在下面提出了 GPR 的六种近似方案，即回归子集 (SR)、Nystrom 方法、数据点子集 (SD)、投影过程 (PP) 近似、贝叶斯委员会机 (BCM) 和迭代解线性系统。 8.3.7 节对这些方法进行了总结，并比较了它们在 2.5 节介绍的 SARCOS 数据上的性能。

8.3.1 回归子集

西尔弗曼 [1985，秒。 6.1]表明可以从有限维广义线性回归模型 f (x*) = ∑n i=1 αik(x*, xi) 中获得均值高斯过程预测因子，其中先验 α ∼ N (0, K−1) .为了看到这一点，我们使用等式给出的特征空间中线性回归模型的平均预测。 (2.11)，即 ̄ f (x∗) = \sigma^{-2}_n φ(x∗)^{\top}A−1Φy 其中 A = Σ−1 p + \sigma^{-2}_n ΦΦ^{\top}。设 φ(x∗) = k(x∗), Φ = Φ^{\top} = K 和 Σ−1 p = K 我们得到 ̄ f (x∗) = \sigma^{-2}_n k^{\top}(x∗)[\sigma^{-2}_n K (K + \sigma^2_nI)]−1Ky (8.11) = k^{\top}(x∗)(K + \sigma^2_nI)^{-1}y, (8.12)

与等式一致。 (2.25)。但是请注意，此模型的预测（协）方差不同于完整的 GPR。

该模型的一个简单近似是仅考虑回归量的一个子集，因此 fSR(x∗) = m ∑ i=1 αik(x∗, xi)，其中 αm ∼ N (0, K^{-1}_{mm})。 (8.13)

再次使用等式。 (2.11) 我们得到̄ fSR(x∗) = km(x∗)^{\top}(K_{mn}K_{nm} + \sigma^2_nK_{mm})^{-1}K_{mn}y, (8.14) V[fSR(x∗)] = \sigma^2_nkm(x∗)^{\top}(K_{mn}K_{nm} + \sigma^2_nK_{mm})^{-1}km(x*)。 (8.15) 因此，αm 的后验均值由 ̄ αm = (K_{mn}K_{nm} + \sigma^2_nK_{mm})^{-1}K_{mn}y 给出。 (8.16)

例如，Wahba [1990，第 7 章] 以及 Poggio 和 Girosi [1990，eq. 25] 通过正则化框架。 G. Wahba 向我们建议了“回归子集”(SR) 这个名称。公式 8.14 和 8.15 的计算需要时间 \mamthcal{O}(m2n) 来执行必要的矩阵计算。此后，新测试点的均值预测需要时间 \mamthcal{O}(m)，预测方差需要 \mamthcal{O}(m2)。

定义为等式。 (8.4)。因此，该模型下的对数边缘似然为 log pSR(y|X) = − 1 2 log | \tilde{K} + \sigma^2_n I_n| − 1 2 y^{\top}( \tilde{K} + \sigma^2_n I_n)^{-1}y − n 2 log(2π)。 (8.17)

请注意，SR 模型定义的协方差函数的形式为 ̃ k(x, x′) = k(x)^{\top}K^{-1}_{mm}k(x′)，这与 Nystrom 近似的形式完全相同协方差函数方程(8.7)。事实上，如果将预测均值和方差方程 2.25 和 2.26 中的协方差函数 k(\mathbf{x,x’}) 系统地替换为 k(\mathbf{x,x’})，我们将得到方程 8.14 和 8.15，如附录 8.6 所示。

如果核函数对于 |x| 衰减为零→ ∞ 对于固定的 x’，那么当 x 远离集合 I 中的点时，̃ k(x, x) 将接近零。即使核是静止的，因此 k(x, x) 是独立于 x。因此，当 x 远离集合 I 中的点时，我们可能期望使用近似核会给出较差的预测，尤其是低估预测方差。

Rasmussen 和 Qui nonero-Candela [2005] 提出的一个缓解这个问题的有趣想法是定义具有 m + 1 个基函数的 SR 模型，其中额外的基函数以测试点 x∗ 为中心，因此 ySR∗ (x*) = ∑m i=1 αik(x*, xi) + αk(x, x*)。然后可以使用该模型进行预测，并且可以使用分区矩阵逆方程 A.11 和 A.12 有效地实现它。以 x* 为中心的额外基函数的作用是保持测试点的预测方差。

到目前为止，我们还没有说应该如何选择子集 I。一种简单的方法是从 X 中随机选择它，另一种方法是在 {xi}n i=1 上运行聚类以获得中心。或者，已经提出了许多 I 的贪婪前向选择算法。 Luo 和 Wahba [1997] 选择下一个核以在优化 αm 后最小化残差平方和 (RSS) |y − K_{nm}αm|2。 Smola 和 Bartlett [2001] 采用了类似的方法，但选择二次型作为他们的标准

\sigma^2_n |y − K_{nm} ̄ αm|2 + ̄ α^{\top} mK_{mm} ̄ αm = y^{\top}( \tilde{K} + \sigma^2_n I_n)^{-1}y, (8.18)

其中右侧跟随使用等式。 (8.16) 和矩阵求逆引理。或者，Qui ̃ nonero-Candela [2004] 建议使用近似对数边缘似然 log pSR(y|X)（参见等式 (8.17)）作为选择标准。事实上，方程式的二次项。 (8.18) 是组成 log pSR(y|X) 的项之一。

对于所有这些建议，通过使用分区矩阵方程，在新示例上评估标准的复杂性为 \mamthcal{O}(mn)。因此，在每次迭代中考虑 R 中的所有点可能过于昂贵，我们可能希望考虑较小的工作集，如算法 8.1 中所述。

请注意，SR 模型是通过以随机或贪婪的方式选择大小为 m 的数据点的一些子集而获得的。 6.6 节中描述的相关向量机 (RVM) 具有类似的风格，因为它会自动选择（以贪婪的方式）在其扩展中使用哪些数据点。但是，请注意一个重要的区别，即 RVM 在 α 上使用对角线先验，而对于 SR 方法，我们有 αm ∼ N (0, K^{-1}_{mm})。

8.3.2 Nystrom 方法

Williams 和 Seeger [2001] 建议通过在均值和方差预测方程 2.25 和 2.26 中用 \tilde{K} 替换矩阵 K 来近似 GPR 方程，并将此称为近似 GPR 的 Nystrom 方法。请注意，在此提案中，协方差函数 k 并未系统地替换为 ̃ k，仅替换了矩阵 K 的出现。对于 SR 模型，执行必要的矩阵计算的时间复杂度为 \mamthcal{O}(m2n)，然后是测试点的预测均值的 \mamthcal{O}(n) 和预测方差的 \mamthcal{O}(mn)。

威廉姆斯等人的实验证据。 [2002] 建议，对于大 m，SR 和 Nystrom 方法具有相似的性能，但对于小 m，Nystrom 方法可能非常差。此外，k 没有系统地被 ̃ k 替换的事实意味着可能会出现尴尬，例如近似预测方差为负。由于这些原因，我们不推荐 Nystrom 方法优于 SR 方法。然而，当 K 的第 (m + 1) 个特征值 λm+1 远小于 \sigma^2_n 时，Nystrom 方法可能有效。

8.3.3 数据点子集

上述回归子集方法近似于预测分布的形式，尤其是预测均值。全样本高斯过程预测器的另一个简单近似是保留高斯过程预测器，但仅在数据大小为 m 的较小子集上。虽然这显然是对数据的浪费，但如果使用 m 个点获得的预测足够准确以满足我们的需求，那么它还是有意义的。

显然，选择将哪些点纳入活动集 I 是有意义的，通常这是通过贪心算法实现的。但是，必须注意所需的计算量，尤其是在每次迭代时考虑 R 的每个成员时。

劳伦斯等人。 [2003] 建议选择下一个点（或站点）包含在活动集中，即最大化微分熵得分 Δj，H[p(fj)] − H[pnew(fj)]，其中 H[p (fj)] 是站点 j ∈ R 处的高斯熵（它是站点 j 处方差的函数，因为后验是高斯分布，参见等式 (A.20)），并且 H[pnew(fj)]是包含站点 j 的观测值后该站点的熵。设包含前 fj 的后验方差为 vj。由于 p(fj|yI , yj) ∝ p(fj|yI )N (yj|fj, \sigma^2) 我们有 (vnew j )^{-1} = v−1 j + σ−2。利用方差为 v 的高斯熵为 log(2πev)/2 这一事实，我们得到 Δj = 1 2 log(1 + vj/\sigma^2)。 (8.19)

∆j 是 vj 的单调函数，因此可以通过选择具有最大方差的站点来最大化它。劳伦斯等人。 [2003] 将他们的方法称为信息向量机 (IVM)

如果编码简单，在单次迭代中计算 R 中所有位点的方差的复杂度为 \mamthcal{O}(m3 + (n − m)m2)，因为我们需要评估等式。 (2.26) 在每个站点（并且 K_{mm} + \sigma^2_nI 的矩阵求逆可以在 \mamthcal{O}(m3) 中完成一次然后存储）。然而，随着我们逐渐增加矩阵 K_{mm} 和 K_{m(n-m)}，实际上每次包含的成本为 \mamthcal{O}(mn)，导致使用大小为 m 的子集时的总体复杂度为 \mamthcal{O}(m2n)。例如，一旦选择了要包含的站点，矩阵 K_{mm} + \sigma^2_nI 就会通过包含额外的行和列而增长。可以使用等式找到此扩展矩阵的逆矩阵。 (A.12) 尽管在数值上使用 Lawrence 等人中描述的 Cholesky 分解方法会更好。 [2003]。该方案在每个步骤中评估所有 j ∈ R 的 ∆j 以选择包含站点。当 m 很小时，这是有意义的，但随着它变大，从 R. Lawrence 等人的子集中选择候选包含站点是有意义的。 [2003] 将此称为随机贪婪选择方法，并就如何选择子集给出了进一步的想法。

微分熵分数 Δj 不是可用于选址的唯一标准。例如，也可以使用信息增益标准 KL(pnew(fj)||p(fj))（参见 Seeger 等人，2003 年）。这里使用贪心选择启发法类似于主动学习的问题，参见例如麦凯 [1992c]。

8.3.4 投影过程近似

SR 方法具有不吸引人的特征，它基于退化的 GP，方程式中给出的有限维模型。 (8.13)。 SD 方法是一种非退化过程模型，但它只使用了 m 个数据点。投影过程 (PP) 近似也是非退化过程模型，但它可以使用所有 n 个数据点。我们称其为投影过程近似，因为它仅表示 m<nm<n 潜在函数值，但通过将 m 个潜在点投影到 n 个维度来计算涉及所有 n 个数据点的可能性。

基本 GPR 算法的一个问题是似然项要求我们具有 n 个训练点的 f 值。但是，假设我们只明确表示这些值中的 m 个，并将它们表示为 fm。然后 R 中剩余的 f 值表示为 fn−m 具有条件分布 p(fn−m|fm)，其均值由 E[fn−m|fm] = K_{(n-m)m}K−1 给出mmfm.3 假设我们用 N (yn−m|E[fn−m|fm], \sigma^2_nI) 替换 R 中点的真实似然项。还包括集合 I 中点的似然贡献，我们有 q(y|fm) = N (y|K_{nm}K^{-1}_{mm}fm, \sigma^2_nI), (8.20)

也可以写成 q(y|fm) = N (y|E[f |fm], \sigma^2_nI)。这里的关键特征是我们将 D 的所有 n 个点中的信息吸收到 I 中的 m 个点中。

方程式中 q(y|fm) 的形式。 (8.20) 可能看起来相当武断，但实际上可以证明，如果我们考虑最小化 KL(q(f |y)||p(f |y))，则近似分布 q(f |y) 之间的 KLdivergence以及形式为 q(f |y) ∝ p(f )R(fm) 的所有 q 分布的真实后验 p(f |y)，其中 R 为正且仅取决于 fm，这就是我们获得的形式。参见 Seeger [2003, Lemma 4.1 and sec. C.2.1] 的详细推导，以及 Csat ´o [2002，秒。 3.3]。

为了进行预测，我们首先必须计算后验分布 q(fm|y)。定义简写 P = K^{-1}{mm}K{mn} 使得 E[f |fm] = P >fm。那么我们有 q(y|fm) ∝ exp ( − 1 2\sigma^2_n (y − P >fm)>(y − P >fm))。 (8.21)

将其与先验 p(fm) ∝ exp(−f > m K −1 mmfm/2) 结合，我们得到 q(fm|y) ∝ exp ( − 1 2f> m (K −1 mm + 1 \sigma^2_n P P > )fm + 1 \sigma^2_n y>P >fm ), (8.22)

可以将其识别为高斯 N (μ, A)，其中 A−1 = \sigma^{-2}n (\sigma^2_n K −1 mm + P P >) = \sigma^{-2}n K^{-1}{mm}(\sigma^2_nK{mm} + K_{mn}K_{nm})K− 1 mm，(8.23) μ = \sigma^{-2}n AP y = K{mm}(\sigma^2_nK_{mm} + K_{mn}K_{nm})^{-1}K_{mn}y。 (8.24)

因此，预测平均值由 Eq[f (x∗)] = km(x∗)>K^{-1}{mm}μ (8.25) = km(x∗)>(\sigma^2_nK{mm} + K_{mn}K_{nm})^{-1}K_{mn}y 给出， (8.26)

事实证明这与 SR 模型下的预测均值相同，如等式 1 所示。 (8.14)。然而，预测方差是不同的。该参数与等式中的相同。 (3.23) 并得出 Vq[f (x∗)] = k(x∗, x∗) − km(x∗)>K^{-1}{mm}km(x∗) + km(x∗)>K^{-1}{mm} cov( fm |y)K −1 mmkm(x∗) = k(x∗, x∗) − km(x∗)>K^{-1}{mm}km(x∗) + \sigma^2_nkm(x∗)>(\sigma^2_nK{mm} + K_{mn}K_{nm}) −1km(x*)。 (8.27)

请注意，预测方差是 SR 模型下的预测方差（等式 (8.27) 中的最后一项）加上 k(x∗, x∗) − km(x∗)>K^{-1}_{mm}km(x∗) 的总和，其中是给定 fm 时 x* 处的预测方差。因此当量。 (8.27) 永远不会小于 SR 预测方差，并且当 x∗ 远离集合 I 中的点时将变得接近 k(x∗, x∗)。

至于 SR 模型，它需要时间 \mamthcal{O}(m2n) 来执行必要的矩阵计算。此后，新测试点的均值预测需要时间 \mamthcal{O}(m)，预测方差需要 \mamthcal{O}(m2)。我们有 q(y|fm) = N (y|P >fm, \sigma^2_nI) 和 p(fm) = N (0, K_{mm})。通过对 fm 进行积分，我们发现 y ∼ N (0, \tilde{K} + \sigma^2_n I_n)。因此，投影过程近似的边缘似然与 SR 模型等式的边缘似然相同。 (8.17)。

再次出现如何选择哪些点进入集合 I 的问题。 Csat o 和 Opper [2002] 提出了一种方法，其中训练示例按顺序呈现（以“在线”方式）。给定当前活动集 I，可以计算新输入点的新颖性；如果这很大，则将此点添加到 I，否则将此点添加到 R。准确地说，输入 x 的新颖性计算为 k(x, x) − km(x)>K^{-1}_{mm}k (x)，在给定 I 中的点的非噪声观察的情况下，可以将其识别为 x 处的预测方差。如果活动集大于某个预设的最大大小，则可以从 I 中删除点，如第 3.3 节中所述Csat ́o 和 Opper [2002]。 Csat o 等人的后期工作。 [2002] 用期望传播类型算法（参见第 3.6 节）取代了上述算法对输入序列的依赖性。

作为选择活动集的替代方法，Seeger 等人。 [2003] 建议按照算法 8.1 使用贪心子集选择方法。合并新站点后信息增益标准的计算需要 \mamthcal{O}(mn)，因此用作选择标准的成本太高。然而，可以廉价地计算信息增益的近似值（请参阅 Seeger 等人 [2003, eq. 3] 和 Seeger [2003, sec. C.4.2] 了解更多详细信息），这允许运行贪心子集算法在每次迭代中 R 中的所有点上。

8.3.5 贝叶斯委员会机

Tresp [2000] 引入了贝叶斯委员会机器 (BCM) 作为加速高斯过程回归的一种方式。令 f* 为测试位置的函数值向量。在 GPR 下，我们获得了 p(f*|D) 的预测高斯分布。对于 BCM，我们将数据集分成 p 个部分 D1, . . . , Dp 其中 Di = (Xi, yi) 并近似 p(y1, . . . , yp|f∗, X) ’ ∠p i=1 p(yi|f∗, Xi)。在这种近似下，我们有 q(f∗|D1, . . . , Dp) ∝ p(f∗) p ‖ i=1 p(yi|f∗, Xi) = c ‖p i=1 p(f∗|Di ) pp−1(f*) , (8.28)

其中 c 是归一化常数。使用分子和面额中的项都是 f* 上的高斯分布这一事实，很容易证明（见练习 8.7.1）f* 的预测均值和协方差由 Eq[f*|D] = [covq(f∗|D)] p ∑ i=1 [cov(f∗|Di)]−1E[f∗|Di], (8.29) [covq(f∗|D)]−1 = −(p − 1)K−1 ∗∗ + p ∑ i=1 [cov(f∗|Di)]−1, (8.30)

其中 K** 是在测试点评估的协方差矩阵。这里 E[f*|Di] 和 cov(f*|Di) 是给定 Di 的 f* 预测的均值和协方差，如等式中给出。 (2.23) 和 (2.24)。请注意，当量。 (8.29) 有一个有趣的形式，因为数据集的每个部分的预测都由逆预测协方差加权。

我们可以自由选择如何对数据集 D 进行分区。这有两个方面，分区的数量和分区的数据点分配。如果我们希望每个分区的大小为 m，则 p = n/m。 Tresp [2000] 使用随机分配数据点到分区，但 Schwaighofer 和 Tresp [2003] 建议对数据进行聚类（例如使用 p-means 聚类）可以提高性能。但是，请注意，与上面使用的贪婪方案相比，聚类不使用目标 y 值，仅使用输入 x。

尽管可以对任意数量的测试点 n* 进行预测，但这会减慢方法的速度，因为它涉及 n* × n* 矩阵的求逆。 Schwaighofer 和 Tresp [2003] 建议对大小为 m 的块进行测试预测，以便所有矩阵的大小相同。在这种情况下，BCM 的计算复杂度是 \mamthcal{O}(pm3) = \mamthcal{O}(m2n) 用于预测 m 个测试点，或每个测试点 \mamthcal{O}(mn)。

BCM 方法是转导的 [Vapnik，1995] 而不是归纳的，因为该方法使用测试集输入位置计算测试集相关模型。另请注意，如果我们希望仅在一个测试点做出预测，则有必要“幻觉”一些额外的测试点作为等式。随着测试点数量的增加，(8.28) 通常会变成更好的近似值。

8.3.6 线性系统的迭代求解

加速高斯过程回归的一种直接方法是注意线性系统 (K + \sigma^2_nI)v = y 可以通过迭代方法求解，例如共轭梯度 (CG)。 （有关 CG 方法的更多详细信息，请参阅 Golub 和 Van Loan [1989，第 10.2 节]。）如果运行 n 次迭代，共轭梯度会给出精确解（忽略舍入误差），但如果终止，它将给出近似解早些时候，假设在 k 次迭代之后，时间复杂度为 \mamthcal{O}(kn2)。 Wahba 等人提出了这种方法。 [1995]（在数值天气预报的背景下）和 Gibbs 和 MacKay [1997]（在一般高斯过程回归的背景下）。 CG 方法也被用于拉普拉斯 GPC 的上下文中，其中重复求解线性系统以获得 MAP 解 f ̃（详见 3.4 和 3.5 节）。

可以加速 CG 方法的一种方法是使用近似而不是精确的矩阵向量乘法。例如，Yang 等人最近的工作。 [2005] 为此目的使用改进的快速高斯变换。

8.3.7 近似 GPR 方法的比较

上面我们介绍了 GPR 的六种近似方法。其中，我们只保留那些与 n 成线性比例关系的方法，因此必须忽略线性系统的迭代解决方案。我们还优先考虑 Nystro ̈m 近似而不是 SR 方法，留下四个备选方案：回归子集 (SR)、数据子集 (SD)、投影过程 (PP) 和贝叶斯委员会机器 (BCM)

表 8.1 显示了 2.5 节中描述的机械臂逆动力学问题的四种方法的结果，其中 D = 21 个输入变量，44,484 个训练示例和 4,449 个测试示例。与第 2.5 节一样，我们对 21 个输入维度中的每一个维度都使用了带有单独长度尺度参数的平方指数协方差函数

对于 SD 方法，随机选择大小为 m 的训练数据子集，并通过优化该子集的边缘似然来设置超参数。由于使用了 ARD，这涉及到 D + 2 个超参数的优化。此过程重复 10 次，产生表 8.1 中记录的平均值和标准偏差。对于 SR、PP 和 BCM 方法，使用了与从 SD 实验中获得的相同的数据和超参数向量子集。4 请注意，m = 4096 结果不适用于 BCM，因为这给出了一个错误的结果内存错误。

这些实验是在具有 3.74 GB RAM 的 2.0 GHz 双处理器机器上进行的。所有四种方法的代码都是用 Matlab 编写的。

表 8.2 总结了四种方法的时间复杂度。因此，对于大小为 n* 的测试集并使用完整（均值和方差）预测，我们发现 SD 方法的时间复杂度为 \mamthcal{O}(m3) + \mamthcal{O}(m2n*)，对于 SR 和 PP 方法，时间复杂度为 \mamthcal{O}(m2n) + \mamthcal{O}(m2n*)，对于 BCM 方法，它是 \mamthcal{O}(mnn*)。假设 n∗ ≥ m 这些分别减少到 \mamthcal{O}(m2n∗)、\mamthcal{O}(m2n) 和 \mamthcal{O}(mnn∗)。这些复杂性与表 8.1 中的时间安排大体一致。

表 8.1 的结果绘制在图 8.1 中。正如我们所料，总的趋势是随着 m 的增加，SMSE 和 MSLL 分数会降低。请注意，使用较小的 m 运行以获得关于 m 的学习曲线是非常值得的；这有助于了解大 m 运行的有用性。就 SMSE 而言，我们看到（毫不奇怪）SD 不如其他方法，而其他方法都具有相似的性能。对于 MSLL，SD 再次不如其他方法，尽管此处 PP 方法对于较大的 m 不如 SR 和 BCM。

这些结果是使用随机选择的活动集获得的。还对 SD 方法 (IVM) 和 SR 方法使用主动选择进行了一些实验，但这些实验并没有显着提高性能。对于 BCM，我们还尝试使用 p 均值聚类而不是随机分配给分区；同样，这并没有导致性能的显着提高。总的来说，在这个数据集上，我们的结论是，对于固定的 m，SR 是首选方法，因为 BCM 对于相似的性能具有更长的运行时间。但是，请注意，如果我们在运行时进行比较，那么 m = 4096 的 SD 在时间和性能上都与 m = 1024 的 SR 和 BCM 结果具有竞争力。

在上面的实验中，所有方法的超参数都是通过优化大小为 m 的 SD 模型的边缘似然来设置的。这意味着我们可以直接比较使用相同超参数和子集的不同方法。但是，可以选择优化每种方法的（近似）边缘似然（参见第 8.5 节），然后比较结果。请注意，优化近似边缘似然的超参数可能取决于方法。例如，图 5.3(b) 显示随着数据量的增加，边缘似然的最大值出现在较短的长度尺度上。 V. Tresp 和 A. Schwaighofer（个人通讯，2004 年）在比较 SD 边缘似然方程时也观察到了这种效应。 (8.31) 对所有 n 个数据点等式计算完全边缘似然。 (5.8)。

Schwaighofer 和 Tresp [2003] 报告了 BCM 方法与一些其他近似方法之间针对许多综合回归问题的一些实验比较。在这些实验中，他们分别为每种方法优化了核超参数。他们的结果是，对于固定的 m，BCM 的性能与其他方法一样好或更好。然而，这些结果取决于数据生成过程中的噪声水平等因素；他们报告（pers. comm., 2005）对于相对较大的噪音水平，BCM 不再显示出优势。根据目前可用的证据，我们无法就一种近似方法优于另一种方法提供坚定的建议；需要进一步研究以了解影响性能的因素。

8.4 具有固定超参数的 GPC 的近似

GPC 的近似方法与 GPR 的近似方法类似，但也需要处理非高斯似然，可以通过使用拉普拉斯近似（请参见第 3.4 节）或期望传播 (EP)（请参见第 3.6 节）。在本节中，我们主要关注二元分类任务，尽管一些方法也可以扩展到多类情况。

对于回归子集 (SR) 方法，我们再次使用模型 fSR(x*) = ∑m i=1 αik(x*, xi) 和 αm ∼ N (0, K^{-1}_{mm})。可能性是非高斯分布的，但寻找 αm 的 MAP 值的优化问题是凸的，可以使用牛顿迭代获得。在这一点上使用 MAP 值 ^ αm 和 Hessian 矩阵，我们获得了 f (x*) 的预测均值和方差，可以通过 sigmoid 函数将其馈入以产生概率预测。像往常一样，出现了如何选择点子集的问题；林等。 [2000] 使用聚类方法选择这些，而 Zhu 和 Hastie [2002] 提出了前向选择策略。

Lawrence 等人提出了用于 GPC 的数据点子集 (SD) 方法。 [2003]，使用后验的 EP 式近似和微分熵得分（参见第 8.3.3 节）来选择要包含的新站点。请注意，EP 近似非常自然地适用于稀疏化：当某些位置精度（参见等式（3.51））为零时会产生稀疏模型，从而使相应的似然项消失。因此，可以通过忽略位置精度非常小的似然项来实现计算增益。

投影过程 (PP) 近似也可以用于非高斯似然。 Csat o 和 Opper [2002] 提出了一种“在线”方法，其中示例按顺序处理，而 Csat o 等人。 [2002] 给出了一种期望传播类型的算法，其中允许多次扫描训练数据。

贝叶斯委员会机器 (BCM) 也被推广到处理 Tresp [2000] 中的非高斯可能性。与 GPR 情况一样，数据集被分成块，但现在在每个块中使用拉普拉斯近似进行近似推理，以产生近似预测均值 Eq[f∗|Di] 和近似预测协方差 covq(f∗|Di ).然后像以前一样使用方程 8.29 和 8.30 组合这些预测。

8.5 近似边缘似然及其导数

我们首先考虑高斯过程回归的近似值，然后是高斯过程分类。对于 GPR，SR 和 PP 方法都会产生与方程式中给出的相同的近似边缘似然。 (8.17)。对于 SD 方法，一个非常简单的近似值（忽略不在活动集中的数据点）由 log pSD(ym|Xm) = − 1 2 log |K_{mm} + \sigma^2I| 给出。 − 1 2 y> m(K_{mm} + \sigma^2I)^{-1}ym − m 2 log(2π), (8.31)

其中 ym 是 y 对应于活动集的子向量；当量。 (8.31) 只是模型 ym ∼ N (0, K_{mm} + \sigma^2I) 下的对数边缘似然。

对于 BCM，一个简单的方法是求和 eq。 (8.31) 在数据集的每个分区上进行评估。这忽略了分区之间的交互。 Tresp 和 Schwaighofer (pers. comm., 2004) 提出了一种更复杂的基于 BCM 的方法，该方法大致考虑了这些相互作用。

对于 SR 近似下的 GPC，可以简单地在有限维模型上使用拉普拉斯或 EP 近似。对于 SD，可以再次忽略所有不在活动集中的数据点，并使用拉普拉斯或 EP 计算对数 p(ym|Xm) 的近似值。对于计划过程 (PP) 方法，Seeger [2003，p. 162] 建议以下下限 log p(y|X) = log ∫ p(y|f )p(f ) df = log ∫ q(f ) p(y|f )p(f ) q(f ) df ≥ ∫ q(f ) log ( p(y|f )p(f ) q(f ) ) df (8.32) = ∫ q(f ) log q(y|f ) df − KL(q(f )|| p(f )) = n ∑ i=1 ∫ q(fi) log p(yi|fi) dfi − KL(q(fm)||p(fm)),

其中 q(f ) 是 q(f |y) 和 eq 的简写。 (8.32) 使用 Jensen 不等式从上一行的等式推导出来。 KL 散度项可以很容易地使用等式进行评估。 (A.23)，并且一维积分可以使用数值求积来解决。

我们不知道将 BCM 近似值扩展到 GPC 的边缘可能性的工作。

鉴于上述边缘似然的各种近似，我们可能还想计算导数以优化它。显然，在优化期间保持活动集固定是有意义的，但请注意，这与选择活动集的方法可能会随着协方差函数参数 θ 的变化而选择不同的集这一事实相冲突。对于分类情况，导数可能非常复杂，因为站点参数（例如 MAP 值 f^，请参阅第 3.4.1 节）随着 θ 的变化而变化。 （我们已经在第 5.5 节中看到了关于非稀疏拉普拉斯近似的示例。）Seeger [2003，sec. 4.8] 描述了一些比较 SD 和 PP 方法优化回归和分类问题边缘似然的实验。

8.6 附录：使用 Nystrom ̈om 近似核的 SR 和 GPR 等价

在第 8.3 节中，我们推导出均值和方差的回归预测子集，如等式 8.14 和 8.15 中给出。本附录的目的是表明这些等同于通过在 GPR 预测方程 2.25 和 2.26 中系统地用 ̃ k(\mathbf{x,x’}) 替换 k(\mathbf{x,x’}) 获得的预测变量。

先为均值。 GPR 预测器是 E[f (x∗)] = k(x∗)>(K + \sigma^2_n I )^{-1} y。用 ̃ k(x, x′) 替换所有出现的 k(x, x′) 我们得到 E[ ̃ f (x∗)] = ̃ k(x∗)>( \tilde{K} + \sigma^2_nI)^{-1}y (8.33 ) = km(x∗)>K^{-1}{mm} K{mn} (K_{nm} K −1 mmK_{mn} + \sigma^2_nI)^{-1}y (8.34) = \sigma^{-2}n km(x∗)>K^{-1}{mm}K_{mn} [ I_n − K_{nm}Q−1K_{mn} ] y (8.35) = \sigma^{-2}n km(x∗)>K^{-1}{mm} [Im − K_{mn}K_{nm}Q−1] K_{mn}y (8.36) = \sigma^{-2}n km(x∗)>K^{-1}{mm} [\sigma^2_nK_{mm}Q−1 ] K_{mn}y (8.37) = km(x*)>Q−1K_{mn}y, (8.38)

其中 Q = \sigma^2_nK_{mm} + K_{mn}K_{nm}，与等式一致。 (8.14)。方程式 (8.35) 由方程式得出。 (8.34) 通过使用矩阵求逆引理 eq。 （A.9）和等式。 (8.38) 从等式得出。 (8.36) 使用 Im = (\sigma^2_nK_{mm} + K_{mn}K_{nm})Q−1。对于预测方差，我们有 V[ ̃ f∗] = ̃ k(x∗, x∗) − ̃ k(x∗)>( \tilde{K} + \sigma^2_nI)^{-1} ̃ k(x∗) (8.39) = km (x∗)>K^{-1}{mm}km(x∗)− (8.40) km (x∗ )> K −1 mm K{mn} (K_{nm} K −1 mmK_{mn} + \sigma^2_n I )^{-1} K_{nm} K −1 mmkm(x∗) = km(x*)>K^{-1}{mm}km(x*) − km(x*)>Q−1K{mn}K_{nm}K^{-1}{mm}km(x*) (8.41) = km(x*)> [Im − Q−1K{mn}K_{nm} ] K −1 mmkm(x∗) (8.42) = km(x∗)>Q−1\sigma^2_n K_{mm} K −1 mmkm(x∗) (8.43) = \sigma^2_nkm(x∗)>Q−1km(x∗), ( 8.44)

与等式一致。 (8.15)。方程之间的步骤。 (8.40) 和 (8.41) 是从等式获得的。上面的 (8.34) 和 (8.38)，以及等式。 (8.43) 从等式得出。 (8.42) 使用 Im = (\sigma^2_nK_{mm} + K_{mn}K_{nm})Q−1。