高斯过程（2）：噪声、正则化与核的线性组合

小象实战讲义 · 机器学习算法精讲

上一节我们用「三点加一维」的直观例子推出了高斯过程的预测公式。这一节先把这个模型正式定义出来（原来讲师说「管它叫 f 不叫 y」，谜底就在这里：建模的对象是函数空间里的 $f(x)$ 本身），再把它从无噪声推广到含噪声的观测——你会发现噪声项的角色像一个正则化；最后把一维情形展开看个底朝天：$\mu_*$ 原来就是核函数值的一个线性组合。收尾是三个大问题——打油井、海上搜救、犯罪热点——它们都是「为什么非要一个概率分布、而不只是一个点估计」的最好注脚。

💡 核心导读

• GP 的正式定义：$f(x) \sim \mathcal{GP}(m(x),\ k(x, x'))$——在函数空间里构造的正态分布，均值函数 $m(x)$、协方差函数（核）$k(x, x')$。
• 无噪声到含噪声：$y = f(x) + \epsilon$、$\epsilon \sim N(0, \sigma_y^2)$——协方差矩阵变成 $K_y = K + \sigma_y^2 I$，对角线加上噪声方差，某种程度上就像一个 regularization。
• 一维展开的真面目：$\mu_* = \sum_i \alpha_i\, k(x_i, x^*)$——预测均值是核值的线性组合，权重 $\alpha = K_y^{-1} y$ 由全部数据全局解出。
• 对比 KDE：KDE 是「一个点与所有点的 contribution 线性组合」；GP 是「一个正态分布到另一个维度正态分布之间的边缘概率关系」。
• 无参的立场：希望数据对自己的行为负责，尽量减少人给定的先验模型对数据本身的影响。

把 GP 正式定义出来

回过头来，把上一节那个「三点加一维」的预测问题重新 refine 一下：所谓 GP，就是我们认为 $x$ 是随机变量，它有一个可构造的函数空间 $f(x)$，而 $f(x)$ 服从一个 Gaussian process——

$$f(x) \sim \mathcal{GP}\big(m(x),\ k(x, x')\big)$$

它以 $m(x)$ 为中心、以 $k(x, x')$ 为协方差。$m(x)$ 是关于 $f(x)$ 的均值函数；$k(x, x')$ 是协方差函数——按定义，协方差是每一个维度的变量都减去自己的均值，再乘到一起，取期望：

$$k(x, x') = \mathbb{E}\big[(f(x) - m(x))(f(x') - m(x'))\big]$$

（方差就是每个点减去均值自己乘自己；协方差是两个不同变量各自减去自己的均值再相乘——这正是 2.3 节学过的定义在函数空间的版本。）

这里要明确一件事：我们一般说的「随机变量服从正态分布」，是变量的取值服从一个分布，$x \sim N(\mu, \sigma)$；而高斯过程，是在函数空间里构造的一个正态分布——是 $x$ 所对应的函数，在函数空间里服从正态分布。这就是 GP 跟别的模型很不一样的出发点。

泛化：noise-free 模型与 N* 个测试点

刚才只是取了个极端的例子（给了三个点、加一个点），现在把它泛化成更一般的形式。第一个先说 noise-free model（无噪声模型）：给定的训练数据 $\lbrace (x_i, f_i)\rbrace_{i=1}^n$，假设 $x_i$ 和 $f_i$ 之间没有任何噪音，满足一种完整的函数关系（它可能通过某种函数关系构造，但是什么我们并不关心），只知道它在函数空间里服从正态分布。

上一节的 $x^*$ 只属于一维的 $\mathbb{R}$；现在推广：$x^*$ 可以属于 $\mathbb{R}^{N_*}$——一次放三个、四个测试点都行，$N_*$ 就是测试点的个数。把联合分布的公式重新写：

$$\begin{bmatrix} \bar{f} \\ f(x^*) \end{bmatrix} \sim \mathcal{N}\left(\begin{bmatrix} m(X) \\ m(x^*) \end{bmatrix},\ \begin{bmatrix} \mathbf{K} & \mathbf{k}_* \\ \mathbf{k}_*^\top & \mathbf{K}_{**} \end{bmatrix}\right)$$

维度的账要算清楚：$\mathbf{K}$ 是 $n \times n$；$\mathbf{k}_*$ 是 $n \times N_*$（$N^*=1$ 时只加一列一行，$N^*=2$ 加两列两行，$N^*=3$ 加三列三行）；$\mathbf{K}_{**}$ 是 $N_* \times N_*$。而其中任何一个 $k$ 里的关系，都满足 kernel function——本节的写法带一个幅度系数：

$$k(a, b) = \sigma_f^2 \exp\left(-\frac{(a - b)^2}{2\ell^2}\right)$$

（这里有两个各自独立的量：$\sigma_f^2$ 是输出幅度的方差系数，决定核在对角线上的取值；$\ell$ 是长度尺度（length scale），控制两点相距多远之后相关性衰减。上一节的 $\lambda$ 写法对应的是长度尺度那一半：$\lambda = 1/(2\ell^2)$。讲师现场把幅度与长度尺度用了同一个符号，本节按标准记法区分开。）

有了这步就可以求解了——每一步得到的结果，都是先验概率到后验概率的推导：想求 $p(f_*)$，给定 $x^*$、再给定训练数据 $X$ 和 $f$，推出来的是一个正态分布——以 $\mu_*$ 为基础、$\sigma_*$ 为方差（或方差矩阵）的正态分布。无论 $f_*$ 有多大的维度：一维时 $\mu_*$ 和 $\sigma_*$ 都是数；多维时 $\mu_*$ 是向量、$\sigma_*$ 是矩阵（这里不详细展开，理解这件事就行）。最终公式，也就是上一节结果加上均值项的泛化版：

$$\mu_* = m(x^*) + \mathbf{k}_*^\top \mathbf{K}^{-1} (\bar{f} - m(X)), \qquad \sigma_*^2 = K_{**} - \mathbf{k}_*^\top \mathbf{K}^{-1} \mathbf{k}_*$$

跟刚才描述的是一样的，只是写成了更泛化的形式。做高斯过程，本质就是假设函数空间在高维空间里呈正态分布，里边的 $x$ 对应的随机变量就是一个 GP，所有重要关系都通过协方差矩阵定义——于是可以根据 $X$ 之间的关系，估计出任何一个 $x^*$ 所对应的概率密度。

GP 与 KDE：边缘概率 vs 线性组合

这里讲师把两节课的差别点破了。KDE 的时候，我们想象的是：其中一个点，是所有其他跟它有关系的点、考虑到它们 contribution 的一个线性组合。而 GP 考虑的是：一个正态分布，到另一个维度的正态分布之间的边缘概率关系——高维空间中呈一个正态分布，它的子空间也是一个正态分布，于是可以推出所有边缘概率之间的关系，通过最开始的 $\mathbf{K}$ 的空间和加上去的 $\mathbf{k}_*$ 的空间，构造出 $f_*$ 的概率密度——这个概率密度，就是对未知变量的一个很好的估计。

一个是线性组合，一个是边缘概率的推导——这是两种思路真正的分野，值得仔细想一下。（不过先别急着下结论，本节后面把一维情形展开时还会再回到这里。）

为什么要一个分布：油井、搜救与犯罪

为什么需要学一个概率分布，而不只是一个点的估计？很多情况下，只对一个点做估计或分类确实够用——比如 density classification：给定一个点，估计它属于 $C_1$ 还是 $C_2$ 的概率谁大谁小，就可以分类判断了。但对连续变量的 prediction，会发现我们除了给出最大估计之外，还给得了相应的概率分布——而这个分布非常有用。

打油井的例子。 比如在新疆或哪个油田（讲师在这里换了块新黑板）：二维平面上，这里打出来一口油井、那里也打出来一口、这几处都出了油。注意它对应的关系：我们刚才一直举的都是低维空间的例子，而油井给定的是二维平面上的点，要拟合的是三维空间里的面。可以假设这些油井满足某种分布——什么分布我们并不知道，是一种很神奇的分布。但有了分布之后，想估计平面上任何其他一个点出油的概率多大，就可以画出相应的面——比如一整套等高线。

为什么需要整个分布、而不只是「最高值那一个点」？因为打每一口井是有成本的：地域有分布——这块是平地（plane），那块是山（mountain）或者沙漠。沙漠里打井成本很高，平地上成本低。有了每个点的概率分布，就可以结合不同位置（沙漠/平地）所占的面积算实际成本。于是会出现这样的局面：两个候选点打出油的概率一样大，但一个点近一半的概率落在平原上、另一个落在沙漠里——成本差了好多。只有知道完整的分布，才能做出这种「概率 × 成本」的综合决策。

海上搜救的例子。 飞机在海上失事了，需要在平面上搜索：已知几个关键点都没有搜到——这些「没搜到」本身也是信息，它们把「残骸不在哪里」的概率关系画了出来。剩下的力量该投向哪，才能让搜到残骸的概率增大？这类高维空间的拟合问题，都可以通过 GP 或者刚才说的 KDE 去做、去尝试。

含上噪声：y = f(x) + ε

刚才的关键假设叫 noise free。现在回过头说：如果包含噪声呢？经常的情况是——$x$ 通过完整函数关系映射成 $f$，再加上一些噪音，才变成我们实际观察到的 $y$：

$$y = f(x) + \epsilon, \qquad \epsilon \sim N(0, \sigma_y^2)$$

$\epsilon$ 就是噪声。（这也回答了上一节「为什么管它叫 f 不叫 y」：$f$ 是理想的无噪声函数，$y$ 才是带噪声的观测——现在 notation 回到了 $y$。）实际生活里观察到的数据模型经常就是这样。

有了这个模型，得到 GP 的方式是类似的：估计的部分不再是 $\bar{f}$，而是 $\bar{y}$；再加一个新的 $y^*$，服从高维正态（均值 $\mu = 0$ 就不写了）。协方差矩阵的分块几乎全部照旧——$\mathbf{k}_*$、$\mathbf{k}_*^\top$、$K_{**}$ 都不变，唯一变的是左上角：

$$\mathbf{K}_y = \mathbf{K} + \sigma_y^2 \mathbf{I}$$

相当于在所有 $k$（$K_{11}, K_{12}, K_{13} \cdots$）的基础上，每一个对角位置上都加上一个噪音。$K_y$ 清楚之后，最终结果把上一节公式里的 $\mathbf{K}$ 全换成 $\mathbf{K}_y$、$\bar{f}$ 换成 $\bar{y}$：

$$\mu_* = \mathbf{k}_*^\top \mathbf{K}_y^{-1} \bar{y}, \qquad \sigma_*^2 = K_{**} - \mathbf{k}_*^\top \mathbf{K}_y^{-1} \mathbf{k}_*$$

把 K 变成 $K_y$，就是多加了一个 $\sigma_y^2$——给定的噪音方差。而你会发现：这个噪音方差在某种程度上，就像一个 regularization 一样（第九章岭回归的解 $(X^\top X + \lambda I)^{-1} X^\top y$ 跟这里 $(K + \sigma_y^2 I)^{-1}$ 是同一个数学动作）。

先看无噪声版本的完整计算，把「线性组合」的真面目一起跑出来（数据为讲义设定，核取 $\sigma_f = 1$、$\ell = 1$）：

import math

def kern(a, b, sf=1.0, ell=1.0):
    """核函数: k(a,b) = sigma_f^2 * exp(-(a-b)^2/(2*ell^2))
    sf = 幅度系数 sigma_f；ell = 长度尺度 length scale"""
    return sf**2 * math.exp(-((a - b) ** 2) / (2 * ell**2))

def inv3(m):
    """3x3 矩阵求逆(伴随矩阵法)"""
    a, b, c = m[0]
    d, e, f = m[1]
    g, h, i = m[2]
    det = a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)
    r = 1.0 / det
    return [[(e*i - f*h)*r, (c*h - b*i)*r, (b*f - c*e)*r],
            [(f*g - d*i)*r, (a*i - c*g)*r, (c*d - a*f)*r],
            [(d*h - e*g)*r, (b*g - a*h)*r, (a*e - b*d)*r]]

def mv(m, v):
    return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))]

X = [0.0, 2.0, 4.0]      # 三个无噪声观测点(讲义设定)
f = [1.0, 4.0, 2.0]
n = 3

K = [[kern(X[i], X[j]) for j in range(n)] for i in range(n)]
alpha = mv(inv3(K), f)   # alpha = K^-1 f
print("无噪声协方差矩阵 K:")
for row in K:
    print("  [" + "  ".join(f"{v:.4f}" for v in row) + "]")
print(f"alpha = K^-1 f = [{', '.join(f'{v:.4f}' for v in alpha)}]")

xq = 1.0
ks = [kern(x, xq) for x in X]
mu = sum(alpha[i] * ks[i] for i in range(n))
print(f"查询点 x* = {xq}, k* = [{', '.join(f'{v:.4f}' for v in ks)}]")
print("mu* 的贡献分解:")
for i in range(n):
    print(f"  alpha[{i}]={alpha[i]:.4f} * k(x[{i}],x*)={ks[i]:.4f} -> {alpha[i]*ks[i]:.4f}")
print(f"mu* = {mu:.4f}")
v_vec = mv(inv3(K), ks)
s2 = kern(xq, xq) - sum(ks[i] * v_vec[i] for i in range(n))
print(f"sigma*^2 = {s2:.4f}")

无噪声协方差矩阵 K:
  [1.0000  0.1353  0.0003]
  [0.1353  1.0000  0.1353]
  [0.0003  0.1353  1.0000]
alpha = K^-1 f = [0.4946, 3.7307, 1.4949]
查询点 x* = 1.0, k* = [0.6065, 0.6065, 0.0111]
mu* 的贡献分解:
  alpha[0]=0.4946 * k(x[0],x*)=0.6065 -> 0.3000
  alpha[1]=3.7307 * k(x[1],x*)=0.6065 -> 2.2628
  alpha[2]=1.4949 * k(x[2],x*)=0.0111 -> 0.0166
mu* = 2.5794
sigma*^2 = 0.3481

看 $\alpha = [0.4946,\ 3.7307,\ 1.4949]$ 这组权重——它就是本节的主角。$\mu_* = 0.4946 \cdot k(x_1, x^*) + 3.7307 \cdot k(x_2, x^*) + 1.4949 \cdot k(x_3, x^*)$：预测均值确实是核函数值的线性组合。注意 $\alpha_2 = 3.73$ 比 $f_2 = 4$ 还小一点点、而 $\alpha_1 = 0.49$ 只有 $f_1 = 1$ 的一半——这组权重不是任何单个数据点的值，而是 $\mathbf{K}^{-1}$ 把三个点之间的相关性「解耦」之后解出来的。

现在把噪声加进去，对照着看（观测值换成带噪声的 $y$，$\sigma_y = 0.5$）：

import math

# 数据与核同上一块(讲义设定), 此块自包含重写
def kern(a, b, sf=1.0, ell=1.0):
    return sf**2 * math.exp(-((a - b) ** 2) / (2 * ell**2))

def inv3(m):
    a, b, c = m[0]
    d, e, f = m[1]
    g, h, i = m[2]
    det = a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)
    r = 1.0 / det
    return [[(e*i - f*h)*r, (c*h - b*i)*r, (b*f - c*e)*r],
            [(f*g - d*i)*r, (a*i - c*g)*r, (c*d - a*f)*r],
            [(d*h - e*g)*r, (b*g - a*h)*r, (a*e - b*d)*r]]

def mv(m, v):
    return [sum(m[i][j] * v[j] for j in range(len(v))) for i in range(len(m))]

X = [0.0, 2.0, 4.0]
y = [1.2, 3.8, 2.1]      # 带噪声的观测值(讲义设定)
n = 3
sy = 0.5                  # 噪声标准差 sigma_y

K = [[kern(X[i], X[j]) for j in range(n)] for i in range(n)]
Ky = [[K[i][j] + (sy**2 if i == j else 0.0) for j in range(n)] for i in range(n)]
alpha_y = mv(inv3(Ky), y)
print(f"含噪声协方差矩阵 K_y = K + sigma_y^2 I (sigma_y={sy}):")
for row in Ky:
    print("  [" + "  ".join(f"{v:.4f}" for v in row) + "]")
print(f"alpha_y = K_y^-1 y = [{', '.join(f'{v:.4f}' for v in alpha_y)}]")

for xq in (1.0, 0.0):
    ks = [kern(x, xq) for x in X]
    mu_n = sum(alpha_y[i] * ks[i] for i in range(n))
    s2_n = kern(xq, xq) - sum(ks[i] * v for i, v in enumerate(mv(inv3(Ky), ks)))
    s2_f = kern(xq, xq) - sum(ks[i] * v for i, v in enumerate(mv(inv3(K), ks)))
    print(f"x*={xq:4.1f}: 含噪声 mu*={mu_n:.4f} sigma*^2={s2_n:.4f} | 无噪声 sigma*^2={s2_f:.4f}")

含噪声协方差矩阵 K_y = K + sigma_y^2 I (sigma_y=0.5):
  [1.2500  0.1353  0.0003]
  [0.1353  1.2500  0.1353]
  [0.0003  0.1353  1.2500]
alpha_y = K_y^-1 y = [0.6543, 2.8204, 1.3745]
x*= 1.0: 含噪声 mu*=2.1227 sigma*^2=0.4670 | 无噪声 sigma*^2=0.3481
x*= 0.0: 含噪声 mu*=1.0364 sigma*^2=0.1994 | 无噪声 sigma*^2=0.0000

三处对照读出「正则化」的味道：第一，$\alpha_y = [0.6543, 2.8204, 1.3745]$ 跟无噪声的 $\alpha$ 完全不同——对角线垫了 $0.25$ 之后，权重被重新分配，$3.73$ 收敛到 $2.82$，不再死盯任何一个观测。第二，$x^* = 0$ 恰好是观测点上：无噪声 GP 的方差是 $0.0000$（精确插值、分毫不让），含噪声后变成 $0.1994$——模型承认「观测本身可能有 $\pm 0.5$ 的抖动」，不再强迫曲线穿过每个点。第三，同一查询点 $x^*=1$ 的方差从 $0.3481$ 升到 $0.4670$——承认噪声，就要为它多付一份不确定性。

一维下的真面目与收官

最后做一遍维度走查（讲师在这段里一板一眼地数过）：$\mathbf{k}_*$ 是 $N \times 1$（$n=3$ 时三行一列）；$\mathbf{k}_*^\top$ 是 $1 \times N$；$\mathbf{K}_y^{-1}$ 是 $N \times N$；$\bar{y}$ 是 $N \times 1$——$1 \times N$ 乘 $N \times N$ 乘 $N \times 1$，最后得到的一定是一个数；$\sigma_*^2$ 那部分乘出来也是一个数。把中间这个 $N \times 1$ 的向量 $\alpha = \mathbf{K}_y^{-1} \bar{y}$ 单独拎出来定义（每个分量 $\alpha_i$ 对应一个训练点），那么：

$$\mu_* = \sum_i \alpha_i \, k(x_i, x^*)$$

也就是说——想求 $x^*$ 对应的 $y^*$，它实际上就是核函数值的一个线性组合。这个情况跟以前描述的 KDE 也好、其他情况也好，其实是类似的：关键还是看每一个点到所有其他点的距离，这个距离用一个 kernel 描述、就是一个权重；GP 在这一点上，也是把这个权重做了一个线性组合、把 kernel 的值做了一个线性组合，赋给新的 $y$ 值。在一维的情况下，这些东西都是相通的。当然 GP 的好处是：可以算高维的情况，同时考虑 $y$ 的噪音的存在——非常广泛，是一个 parameter free（无参数） 的模型（讲师口播说「skill free 或者叫 parameter free」，正确的说法是后者）。

再补两个 GT 的应用场景把这一章收尾。犯罪热点与警力布置（讲师提到这是某年数学建模竞赛的原型）：一个城市（比如北京市各个区）很多犯罪点都标出来了，给定一个新的一小块区域，能不能找到它可能发生犯罪概率最大的分布？警察局完全可以根据已知犯罪点的情况布置不同的警力——考虑到未知点可能发案的概率；案件还可以分类型（简单的抢劫、到杀人、到偷盗），按权重重新算警力分配的比例——这就是 model free 的预测和判断。反过来想，如果硬用参数化模型：认为分布服从高斯？不太对。t 分布？也不太对。什么很奇怪的构造出来的概率？——当你的先验不准确的时候，这个问题就很难，准确率也不好。而用 KDE 或 GP：完全可以根据已知数据点，对整个区域做概率估计，再按概率去算 cost 最小的规划——这是比较好的办法。

所以 Gaussian process 和 KDE 这种无参数模型，今天讲到这里。要记住：它们都是 nonparametric model——机器学习里非常重要的领域，到今天还有很多人在研究。因为我们的立场是：希望数据对自己的行为负责，尽量减少人给定先验模型对数据本身的影响。现实生活中的问题——犯罪也好、找飞机残骸也好——都不是简单的概率或其他参数模型所能描述的，这类问题更希望数据自己起到很重要的作用，这些方法就都是比较好的选择。

📝 动手练一练

1. 「噪声 = 正则化」：为什么把 $\sigma_y^2$ 加到协方差矩阵对角线上，效果「像一个 regularization」？从「模型是否强迫穿过每个观测点」的角度说一说，并指出它与第九章岭回归解 $(X^\top X + \lambda I)^{-1} X^\top y$ 在数学结构上的对应。

2. 两点 GP 手算：观测 $(x_1, y_1) = (0, 1)$、$(x_2, y_2) = (2, 3)$（无噪声），核 $k(a, b) = \exp(-(a-b)^2/2)$，查询点 $x^* = 1$。手算 $\mathbf{K}$、$\alpha = \mathbf{K}^{-1} \bar{y}$、$\mu_*$、$\sigma_*^2$。（提示：$e^{-0.5} \approx 0.6065$，$e^{-2} \approx 0.1353$，$e^{-4} \approx 0.0183$。）

本章小结

这一节把高斯过程收了口。第一，正式定义：$f(x) \sim \mathcal{GP}(m(x), k(x, x'))$——在函数空间里构造的正态分布，均值函数加协方差核；训练与测试联合服从高维正态，$\mathbf{k}_*$ 是 $n \times N_*$、$\mathbf{K}_{**}$ 是 $N_* \times N_*$。第二，每一步都是先验到后验的概率推导——这是 GP 与 KDE 的分野：KDE 是 contribution 的线性组合，GP 是正态分布之间的边缘概率关系。第三，为什么要分布：打油井要算「概率 × 成本」（两个出油概率相同的点，成本可以差一半）、海上搜救要利用「没搜到」的负信息。第四，含噪声模型 $y = f(x) + \epsilon$：$\mathbf{K}_y = \mathbf{K} + \sigma_y^2 \mathbf{I}$，噪声方差的角色像一个 regularization——不再强迫穿过每个观测点。第五，一维真面目：$\mu_* = \sum_i \alpha_i k(x_i, x^*)$，核值的线性组合，权重 $\alpha$ 由全部数据全局解出——一维下与 KDE 相通，而 GP 胜在能算高维、能容噪声、parameter free。收官的立场：希望数据对自己的行为负责，尽量减少人给定先验模型对数据本身的影响。

📋 行动清单

• 口述定义：不借助资料，说出 GP 的正式定义（均值函数 + 协方差核）和「函数空间里的正态分布」这句话的含义。
• 手算走查：独立完成练一练第 2 题的两点 GP，再用正文代码思路验证。
• 跑改代码：把含噪声块的 $\sigma_y$ 改成 0.1 和 2.0 各跑一遍，观察 $\alpha_y$ 的收敛幅度与观测点上方差的变化方向。
• 讲一个例子：用自己的话把「打油井为什么需要整个分布而不是一个点」讲给同学听，要求讲到沙漠与平原的成本差。

—— 小象教研组