一、无监督降维漏掉了什么

先把问题造出来。三个类，每类 100 个点，类心沿纵轴排开，间距 1.2。 但类内横向的标准差是 3.0，纵向只有 0.45。 也就是说：类与类的差别在纵轴，可数据的发散在横轴。 两个散度矩阵可以把这个矛盾量化：

二、两个散度矩阵与广义特征值

判别分析要找一个方向 w，让投影之后「类间散得开、类内挤得紧」。 写成比值就是 Fisher 准则：

J(w) = ( w' * Sb * w ) / ( w' * Sw * w )

分子是投影后的类间散度，分母是投影后的类内散度。 J 大意味着三个类的投影点分得很开、各自又很窄——正是好分类想要的样子。 这个比值对 w 的长度不敏感（上下都是二次型），所以只看方向。 对它求极值，得到的是一个广义特征值问题：

Sb * w = lam * Sw * w     也就是     解 (Sw 逆) * Sb 的特征向量

最大的那个 lam 就是最优的 J 值，对应的 w 就是第一判别方向。 实测两个方向上的 J：

三、降到一维之后

光看准则值还不够，直接比分类准确率更实在。 用最近邻（留一）在这三份数据上分类：

四、判别方向最多有几个

判别分析能给出的方向个数有硬上限：C 减 1，C 是类别数。 原因在 Sb 的构造上——它是 C 个类心相对总均值的偏差做外积再加权求和， 而这 C 个偏差之和为零（因为减的就是加权平均），所以秩最多 C 减 1。 这个数据里 C 等于 3，上限是 2。可实际只拿到 1 个：

五、动手 Demo：判别分析：把标签用上，选出来的方向和 PCA 几乎垂直

这段代码造三个沿纵轴排开的类，但让类内横向散得远大于纵向间距； 然后分别用 PCA 和判别分析各找一个方向，比夹角、比 Fisher 准则值、比降到一维之后的最近邻准确率。 重点看三处。第一段那两个散度矩阵的对角元，一个 2756.39 一个 311.51； 第二段两个方向的夹角 89.5 度；第三段 PCA 一维的 0.3333 和判别分析一维的 0.8567。

import numpy as np


def lda_data(n=300, seed=3, sep=1.2, sx=3.0, sy=0.45):
    """三个类沿纵轴排开，但类内横向的发散远大于纵向的间距。"""
    rg = np.random.default_rng(seed)
    cs = np.array([[0.0, -sep], [0.0, 0.0], [0.0, sep]])
    per = n // 3
    X, y = [], []
    for i, c in enumerate(cs):
        X.append(c + rg.normal(0.0, 1.0, (per, 2)) * np.array([sx, sy]))
        y.append(np.full(per, i))
    return np.vstack(X), np.concatenate(y)


def scatter(X, y):
    """类内散度 Sw 与类间散度 Sb。"""
    d = X.shape[1]
    mu = X.mean(axis=0)
    Sw = np.zeros((d, d))
    Sb = np.zeros((d, d))
    for c in np.unique(y):
        Xc = X[y == c]
        mc = Xc.mean(axis=0)
        D = Xc - mc
        Sw += D.T @ D
        Sb += len(Xc) * np.outer(mc - mu, mc - mu)
    return Sw, Sb


def pca_dir(X, k=2):
    """不用标签：按方差从大到小取方向。"""
    Xc = X - X.mean(axis=0)
    C = Xc.T @ Xc / len(X)
    vals, vecs = np.linalg.eigh(C)
    order = np.argsort(-vals)
    return vecs[:, order][:, :k], vals[order]


def lda_dir(X, y, k=2):
    """用标签：解广义特征值问题 Sb w 等于 lam 乘 Sw w。"""
    Sw, Sb = scatter(X, y)
    vals, vecs = np.linalg.eig(np.linalg.solve(Sw, Sb))
    order = np.argsort(-vals.real)
    return vecs[:, order].real[:, :k], vals[order].real


def fisher_val(X, w, y):
    """方向 w 上的 Fisher 准则：类间散度除以类内散度。"""
    z = X @ w
    mu = z.mean()
    sb, sw = 0.0, 0.0
    for c in np.unique(y):
        zc = z[y == c]
        sb += len(zc) * (zc.mean() - mu) ** 2
        sw += float(((zc - zc.mean()) ** 2).sum())
    return sb / sw if sw > 0 else float("inf")


def nn_acc(Z, y):
    """最近邻分类的留一准确率。"""
    D2 = ((Z[:, None, :] - Z[None, :, :]) ** 2).sum(-1)
    np.fill_diagonal(D2, np.inf)
    return float(np.mean(y[np.argmin(D2, axis=1)] == y))


X, Y = lda_data(300, 3)
Sw, Sb = scatter(X, Y)
Vp, lam_p = pca_dir(X, 2)
Vl, lam_l = lda_dir(X, Y, 2)

print('[1] 数据：三类各 100 个点')
print('  类内散度 Sw = [[%.2f, %.2f], [%.2f, %.2f]]' % (Sw[0, 0], Sw[0, 1], Sw[1, 0], Sw[1, 1]))
print('  类间散度 Sb = [[%.2f, %.2f], [%.2f, %.2f]]' % (Sb[0, 0], Sb[0, 1], Sb[1, 0], Sb[1, 1]))
print('  Sw 的迹 %.2f，Sb 的迹 %.2f——类内的发散是类间间距的 %.1f 倍'
      % (np.trace(Sw), np.trace(Sb), np.trace(Sw) / np.trace(Sb)))

print('')
print('[2] 不用标签和用标签，选出的方向几乎垂直')
print('  PCA 第一方向 (%.4f, %.4f)  方差 %.4f'
      % (Vp[0, 0], Vp[1, 0], lam_p[0]))
print('  LDA 第一方向 (%.4f, %.4f)  广义特征值 %.4f'
      % (Vl[0, 0], Vl[1, 0], lam_l[0]))
c = abs(float(Vp[:, 0] @ Vl[:, 0]))
print('  两个方向的夹角 %.1f 度' % np.degrees(np.arccos(min(c, 1.0))))
print('  Fisher 值：PCA 第一方向 %.4f   LDA 第一方向 %.4f'
      % (fisher_val(X, Vp[:, 0], Y), fisher_val(X, Vl[:, 0], Y)))

print('')
print('[3] 降到一维之后，最近邻分类的留一准确率')
print('  原始二维 %.4f' % nn_acc(X, Y))
print('  取 PCA 第一维 %.4f' % nn_acc(X @ Vp[:, :1], Y))
print('  取 LDA 第一维 %.4f' % nn_acc(X @ Vl[:, :1], Y))

print('')
print('[4] 判别方向最多 C 减 1 个，这里连 C 减 1 都不到')
print('  广义特征值 ' + '  '.join('%.4f' % v for v in lam_l))
ev = np.linalg.eigvalsh(Sb)
print('  Sb 的两个特征值 %.4f 和 %.2e' % (ev[1], ev[0]))
print('  两个特征值差 %.0f 倍：三个类心排在同一条直线上，类间散度实际只有一个方向'
      % (ev[1] / ev[0]))

实跑输出：

--------------------------------------------------------------------
判别分析：把标签用上，选出来的方向和 PCA 几乎垂直
--------------------------------------------------------------------
[1] 数据：三类各 100 个点
  类内散度 Sw = [[2756.39, 12.21], [12.21, 56.56]]
  类间散度 Sb = [[1.13, -17.66], [-17.66, 311.51]]
  Sw 的迹 2812.95，Sb 的迹 312.64——类内的发散是类间间距的 9.0 倍

[2] 不用标签和用标签，选出的方向几乎垂直
  PCA 第一方向 (-1.0000, 0.0023)  方差 9.1918
  LDA 第一方向 (0.0056, -1.0000)  广义特征值 5.5160
  两个方向的夹角 89.5 度
  Fisher 值：PCA 第一方向 0.0004   LDA 第一方向 5.5160

[3] 降到一维之后，最近邻分类的留一准确率
  原始二维 0.8300
  取 PCA 第一维 0.3333
  取 LDA 第一维 0.8567

[4] 判别方向最多 C 减 1 个，这里连 C 减 1 都不到
  广义特征值 5.5160  0.0000
  Sb 的两个特征值 312.5131 和 1.30e-01
  两个特征值差 2403 倍：三个类心排在同一条直线上，类间散度实际只有一个方向

第一段两个矩阵要横着看。 Sw 的左上 2756.39 是横轴方向的类内散度，Sb 的右下 311.51 是纵轴方向的类间散度。 一个管噪声，一个管信号，恰好错开了维度。 Sw 的迹是 Sb 的 9.0 倍，所以只要按总方差排序，横轴一定赢。 第二段两个方向：PCA 给出 (-1.0000, 0.0023)，判别分析给出 (0.0056, -1.0000)。 夹角 89.5 度。 Fisher 准则值 0.0004 对 5.5160——这两个数放在一起， 就是「无监督降维为什么会在分类任务上翻车」的完整证据。 第三段三个准确率 0.8300、0.3333、0.8567。 中间的 0.3333 等于三类的先验概率，说明投影之后最近邻完全失效。 最后的 0.8567 比原始二维还高一点——丢掉噪声方向带来的收益，大过了丢掉信息带来的损失。 第四段两行是在提醒上界：广义特征值 5.5160 和 0.0000。 第二个是零，说明这个数据只有一维可降。 三个类心共线，Sb 的两个特征值差 2403 倍，第二个方向根本没有能量。 代码里有一处容易写错：解广义特征值问题用的是 eig(solve(Sw, Sb))， 而不是 eig(Sb, Sw) 之类。 注意 np.linalg.eig 返回的特征向量不是正交的， 因为这些向量关于 Sw 正交而不是关于欧氏内积正交。 拿去做投影完全没问题，但别拿它们当正交基用。

30 秒实验

1. 把 lda_data 的 sx 从 3.0 改成 0.4，让横向不再那么散，重跑。 这时 PCA 和判别分析会开始接近——两个判据不再打架。 看夹角掉到多少度，就能体会「重要性是相对的」这句话的分量。
2. 把 sx 改成 3.0、sy 改成 1.2，让类内纵向散度和类间间距一样大。 看判别分析一维的准确率掉到多少。 这时候任何一维投影都会重叠，降维本身就变成亏本买卖。
3. 把 scatter 里的 Sw 加一个 1e-3 * np.eye(2) 再解。这是正则判别分析的雏形。 在小样本高维的真实数据上，这个岭项不是可选项，是必需品——不加就解不出来。

常见误区

1. 「PCA 是最优降维，别的都是它的变体」。 PCA 最优的前提是「用方差衡量信息」，这个前提在做分类时并不成立。 实测两个方向夹角 89.5 度，Fisher 准则值差一万多倍。 在分类任务上，PCA 选的方向可以是最差的那个。
2. 「降维一定会掉准确率」。 实测判别分析压到一维是 0.8567，比原始二维的 0.8300 还高。 因为被丢掉的那一维是纯噪声，去掉它等于免费去噪。 降维的收益和损失要一起算。
3. 「判别分析能降到任意维」。 上限是 C 减 1，而且实际可用维数常常更少。 本节 3 个类最多 2 维，实测只有 1 维有能量：Sb 的两个特征值 312.5131 和 0.1300。
4. 「类内协方差不影响判别分析」。 判别分析假设所有类共用一个 Sw，所以它只找得出「平移型」的差别。 如果两个类形状完全不同（比如一个横条一个竖条）， 判别分析会束手无策——那种情况该用二次判别分析，让每个类有自己的协方差。

判别分析有一个前提：它认定「有用的方向」就是「类心分开的方向」。 可如果标签和输入的关系不是靠均值差体现的呢？ 下一节换一个更根本的问法：到底哪几个方向就够描述标签了。

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

相关笔记

• 充分降维
• 主成分分析
• 最小二乘分类
• 🤖 图解机器学习 · 索引