一、把目标函数写成矩阵

先约定记号：全部样本的输入摞成矩阵 X（n 行 d 列），有标签的那一小部分记作 X_l， 对应的标签记作 y_l。 模型是线性的，预测值 f 等于 X w。 目标函数两项：

第一项  || y_l - X_l w || 的平方        # 有标签的要拟合准
第二项  lam * w' * X' * L * X * w       # 相邻点的预测值要接近

为什么第二项长这样？ 因为相邻点的预测差可以写成 X w 这个东西在图上求「差分」， 而求差分的二次型正好是 w 转置乘 X 转置乘 L 乘 X 乘 w。 这一步的推导在 16.3 节会展开，那里会证明它等于「所有边上加权差的平方和的一半」。 对 w 求导：第一项给出 2 乘 X_l 转置乘（X_l w 减 y_l）， 第二项给出 2 lam 乘 X 转置乘 L 乘 X 乘 w。 令导数等于零，移项：

( X_l' X_l + lam * X' L X ) w = X_l' y_l

左边是一个 d 乘 d 的方阵，右边是一个 d 维向量。 解一个线性方程组就完事——用一次矩阵分解，代价和「解 30 个未知数的方程」一样。 对比一下梯度下降要干什么：算 4000 次梯度、每次都要乘一遍 n 乘 d 的矩阵。 闭式解一次就到位，这是拉普拉斯正则化最小二乘最大的工程优势。

二、为什么加了这一项就可解

10 个标签、30 维的时候，X_l 只有 10 行。 X_l 转置乘 X_l 是一个 30 乘 30 的矩阵，但它的秩最多 10——这是一个秩亏矩阵，行列式为零， 不可逆。 实测秩就是 10。 所以 lam 等于 0 时方程有无数个解，程序只能靠加一个极小的岭项强行挑一个， 挑到哪个纯看数值误差。 加上 lam 乘 X 转置乘 L 乘 X 之后，情况变了。 X 转置乘 L 乘 X 是半正定的（L 本身半正定），它把所有样本都用上了，而不只是那 10 个。 只要图连得够好、lam 大于 0，两半加起来就是正定的，解唯一。 换个说法：标签项只约束「10 个点上的取值」，图项约束「所有点之间的相对关系」。 前者给出 10 个约束，后者把剩下的自由度全填满。 这就是半监督在代数层面的真实作用——它补的是秩，不是数据。

三、lam 该取多大

从 0.003 一路扫到 1000，每类 5 个标签、k 等于 8：

四、k 该取多大

固定 lam 等于 1，把 k 从 3 扫到 80：

五、闭式解与迭代解

既然闭式解一次就完，为什么还要谈迭代？ 因为特征维数 d 很大时，解一个 d 乘 d 的方程组会很贵；而梯度下降每次只算一次矩阵乘向量。 这一节的数据 d 等于 30，两个都便宜，正好拿来对比。

六、动手 Demo：拉普拉斯正则化怎么解：一次线性方程组，不用迭代

这段代码先确认「10 个标签时矩阵不满秩」，再扫 lam、扫 k，然后把梯度下降的迭代解和闭式解做对比， 最后看条件数怎么随 lam 变。 重点看三处。第一段的「秩 10」为什么导致 lam 等于 0 时解不出来； 第二段 lam 从 0.003 到 1000 那条先升后平的曲线； 以及第四段四千步还没消掉的 1.9 乘 10 的负 3 次方差距。

import numpy as np


def manifold_data(n=200, d=30, seed=2, noise=0.05):
    """200 个点躺在 30 维空间里的一个二维流形上，标签由线性方向决定。"""
    rg = np.random.default_rng(seed)
    z = rg.normal(0.0, 1.0, (n, 2))
    A = rg.normal(0.0, 1.0, (d, 2))
    X = z @ A.T + rg.normal(0.0, noise, (n, d))
    v = rg.normal(0.0, 1.0, d)
    v = v / np.linalg.norm(v)
    y = np.where(X @ v > 0.0, 1.0, -1.0)
    return X, y, v


def pick_labels(y, per=5, seed=0):
    """每类随机挑 per 个下标当有标签样本。"""
    rg = np.random.default_rng(seed)
    lab = []
    for c in np.unique(y):
        idx = np.where(y == c)[0]
        lab.extend(rg.choice(idx, per, replace=False).tolist())
    return np.array(sorted(lab))


def graph(X, k=8):
    """k 近邻图：W 邻接矩阵，L 等于度矩阵减 W。"""
    n = len(X)
    D2 = np.maximum(((X[:, None, :] - X[None, :, :]) ** 2).sum(-1), 0.0)
    np.fill_diagonal(D2, np.inf)
    idx = np.argsort(D2, axis=1)[:, :k]
    near = np.take_along_axis(D2, idx, axis=1)
    s2 = float(np.median(near))
    W = np.zeros((n, n))
    W[np.repeat(np.arange(n), k), idx.ravel()] = np.exp(-near.ravel() / (2.0 * s2))
    W = np.maximum(W, W.T)
    return W, np.diag(W.sum(axis=1)) - W


def lap_rls(X, y, lab, lam=1.0, k=8):
    """拉普拉斯正则化最小二乘：标签项加上图上相邻点要接近。"""
    W, L = graph(X, k)
    Xl, yl = X[lab], y[lab]
    return np.linalg.solve(Xl.T @ Xl + lam * (X.T @ L @ X), Xl.T @ yl)


def acc(w, X, y):
    return float(np.mean(np.sign(X @ w) == y))


def purity(X, y, k=8):
    """每个点的 k 个近邻里同类的比例，用来看流形假设成不成立。"""
    D2 = np.maximum(((X[:, None, :] - X[None, :, :]) ** 2).sum(-1), 0.0)
    np.fill_diagonal(D2, np.inf)
    idx = np.argsort(D2, axis=1)[:, :k]
    return float((y[idx] == y[:, None]).mean())


X, Y, V = manifold_data(200, 30, 2)
LAB = pick_labels(Y, 5, seed=0)
XL, YL = X[LAB], Y[LAB]

print('[1] 只有 10 个标签，矩阵根本不满秩')
print('  特征维数 30，有标签样本 10 个')
print('  Xl 转置乘 Xl 的秩 %d（满秩需要 30）' % int(np.linalg.matrix_rank(XL.T @ XL)))
print('  所以 lam 等于 0 时解不唯一，必须靠数值兜底才解得出东西')

print('')
print('[2] lam 从小到大（每类 5 个标签，k=8）')
for lm in (0.003, 0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0, 100.0, 1000.0):
    print('  lam=%-9.3f 准确率 %.4f' % (lm, acc(lap_rls(X, Y, LAB, lam=lm, k=8), X, Y)))

print('')
print('[3] k 从小到大（lam=1）')
for k in (3, 5, 8, 12, 20, 40, 80):
    W, L = graph(X, k)
    w = lap_rls(X, Y, LAB, lam=1.0, k=k)
    print('  k=%-3d 准确率 %.4f  有向边 %5d  近邻同类比例 %.4f'
          % (k, acc(w, X, Y), int((W > 0).sum()), purity(X, Y, k)))

print('')
print('[4] 梯度下降能不能走到闭式解（lam=1，k=8）')
W, L = graph(X, 8)
A = XL.T @ XL + 1.0 * (X.T @ L @ X)
B = XL.T @ YL
W_STAR = np.linalg.solve(A, B)
STEP = 1.0 / np.linalg.eigvalsh(A).max()
w = np.zeros(30)
for it in range(1, 4001):
    g = XL.T @ (XL @ w - YL) + X.T @ (L @ (X @ w))
    w = w - STEP * g
    if it in (10, 100, 1000, 2000, 4000):
        print('  第 %4d 步  与闭式解的差距 %.3e  准确率 %.4f'
              % (it, float(np.linalg.norm(w - W_STAR)), acc(w, X, Y)))

print('')
print('[5] 顺带把病态也治了：矩阵条件数')
for lm in (0.0, 0.01, 1.0, 100.0):
    A2 = XL.T @ XL + lm * (X.T @ L @ X) + (0.0 if lm > 0 else 1e-8 * np.eye(30))
    ev = np.linalg.eigvalsh(A2)
    print('  lam=%-6.2f 最小特征值 %.3e  最大 %.3e  条件数 %.3e'
          % (lm, ev.min(), ev.max(), ev.max() / max(ev.min(), 1e-300)))

实跑输出：

--------------------------------------------------------------------
拉普拉斯正则化怎么解：一次线性方程组，不用迭代
--------------------------------------------------------------------
[1] 只有 10 个标签，矩阵根本不满秩
  特征维数 30，有标签样本 10 个
  Xl 转置乘 Xl 的秩 10（满秩需要 30）
  所以 lam 等于 0 时解不唯一，必须靠数值兜底才解得出东西

[2] lam 从小到大（每类 5 个标签，k=8）
  lam=0.003     准确率 0.7950
  lam=0.010     准确率 0.8600
  lam=0.030     准确率 0.9000
  lam=0.100     准确率 0.9550
  lam=0.300     准确率 0.9700
  lam=1.000     准确率 0.9850
  lam=3.000     准确率 1.0000
  lam=10.000    准确率 1.0000
  lam=30.000    准确率 1.0000
  lam=100.000   准确率 1.0000
  lam=1000.000  准确率 0.9950

[3] k 从小到大（lam=1）
  k=3   准确率 0.9550  有向边   754  近邻同类比例 0.9783
  k=5   准确率 0.9850  有向边  1256  近邻同类比例 0.9620
  k=8   准确率 0.9850  有向边  1994  近邻同类比例 0.9500
  k=12  准确率 0.9900  有向边  2992  近邻同类比例 0.9337
  k=20  准确率 0.9750  有向边  4970  近邻同类比例 0.9130
  k=40  准确率 0.9500  有向边  9972  近邻同类比例 0.8614
  k=80  准确率 0.8550  有向边 19872  近邻同类比例 0.7454

[4] 梯度下降能不能走到闭式解（lam=1，k=8）
  第   10 步  与闭式解的差距 2.580e-01  准确率 0.9950
  第  100 步  与闭式解的差距 2.216e-01  准确率 0.9950
  第 1000 步  与闭式解的差距 5.935e-02  准确率 0.9850
  第 2000 步  与闭式解的差距 1.690e-02  准确率 0.9850
  第 4000 步  与闭式解的差距 1.931e-03  准确率 0.9850

[5] 顺带把病态也治了：矩阵条件数
  lam=0.00   最小特征值 1.000e-08  最大 3.404e+02  条件数 3.404e+10
  lam=0.01   最小特征值 1.090e-02  最大 3.484e+02  条件数 3.196e+04
  lam=1.00   最小特征值 8.518e-01  最大 1.147e+03  条件数 1.346e+03
  lam=100.00 最小特征值 8.363e+01  最大 8.166e+04  条件数 9.766e+02

第一段只有三行，但它决定了整节的存在意义。 30 维、10 个标签，X_l 转置乘 X_l 的秩是 10——缺了 20 个维度。 这意味着 lam 等于 0 时方程有无穷多解，程序只能靠 1e-8 的岭项挑一个。 第二段十一行是这条曲线的形状：0.7950、0.8600、0.9000、0.9550、0.9700、0.9850、 然后 3 到 100 全是 1.0000，最后 1000 掉到 0.9950。 从 3 到 100 有三十三倍的区间都是满分，说明这个超参数很好伺候。 第三段七行看两列：准确率先升后降，近邻同类比例一路单调降。 k 等于 12 时 0.9900 最好，k 等于 80 时 0.8550 最差，而最后那一行的近邻同类比例只有 0.7454。 这两个数要一起看：图连错了，正则项就在往错的方向推。 第四段梯度下降的数字有点扎眼：四千步之后差距还有 1.9 乘 10 的负 3 次方。 步长取的是最大特征值的倒数，这是保证不发散的最大步长； 而条件数是 1346，所以收敛速度被压得很慢。 这也说明闭式解的价值：同样的精度，一次分解就够，不用等四千步。 最后一段四行是附带收获：lam 从 0 加到 1，条件数从 3.4e10 降到 1.3e3。 很多人只知道正则化能防过拟合，其实它还有一个身份是「数值稳定器」。

30 秒实验

1. 把 lap_rls 里的 graph(X, k) 换成 graph(X, 1)。 k 等于 1 时每个点只连一个邻居，图基本是断的。 看准确率掉到多少，再想想「图必须连通」这个条件到底有多硬。
2. 把梯度下降的步长改成 2.0 / 最大特征值，重跑第四段。 这是超出稳定上限的步长，应该会震荡甚至发散。 看前几步的差距是变大还是变小。
3. 把 [2] 里的 lam 序列整体乘 10，看那一大段满分平台往哪边挪。 如果平台整体右移，说明最优 lam 和样本量、维数有关，不能照抄。

常见误区

1. 「lam 等于 0 就是用标签的版本，也算一个合法的对照」。 它不是合法的——这个矩阵秩亏，解不唯一。 实测 X_l 转置乘 X_l 的秩只有 10，程序解出来的是「靠岭项随便挑的一个」，换一台机器、 换一个库，结果都可能不同。
2. 「lam 越大越好，因为正则化总能防过拟合」。 lam 到 1000 时准确率从 1.0000 掉到 0.9950。 原因不是过拟合，是欠拟合：平滑项压过标签项之后，模型干脆不去拟合标签了。
3. 「闭式解和梯度下降解出来的是同一个东西，用哪个都行」。 数学上同一个，数值上差很多。 实测梯度下降四千步后与闭式解还差 1.9e-3，条件数 1346 让收敛极慢； 而 d 很大时闭式解又要分解一个 d 乘 d 的矩阵。 选哪个取决于 d 和 n 谁大。
4. 「k 大一点图更连通，信息传得更远，所以更好」。 k 从 12 加到 80，近邻同类比例从 0.9337 掉到 0.7454，准确率从 0.9900 掉到 0.8550。 图连通和信息正确是两件事，跨类的边只会把错的信息传得更远。

解出来了，但那一项到底在惩罚什么，还只是一句「相邻点要接近」。 下一节把它拆开算给你看，顺便给出一个反例：平滑度一路变好，准确率却一路变差。

• ☐ （本节完，下一节 16.3 16.3）

相关笔记

• 灵活应用输入数据的流形构造
• 拉普拉斯正则化的解释
• 最小二乘解的性质
• 🤖 图解机器学习 · 索引