一、线性投影的死角

先把问题摆清楚。400 个二维点，内圈半径 1、外圈半径 3，各 200 个。 要求：只用一个数把两圈分开。 线性 PCA 降到一维，就是找一条穿过圆心的直线，把每个点投影上去。 投影结果是：内圈的点落在中间一段，外圈的点落在两头——因为外圈半径大，投影出来的绝对值也大。 但注意「两头」是分居正负两侧的，中间一段必然和外圈的一部分重叠。 所以任何一刀都切不干净。实测扫遍所有可能阈值，最好的一刀在 0.9078，正确率 72.50%。 这不是阈值没调好，是投影把「半径」这个信息扔掉了。 要切开同心圆，需要的坐标是「半径」——它是原始坐标的平方和开根，不是任何线性组合。 这就是线性方法的死角：它只能输出 x 和 y 的加权和，永远做不出 r 等于根号下 x² 加 y²。 顺带说一句：这里是二维降到一维，所以「直线」就是一条过原点的线。 三维里降到二维，「切不开」的现象是一样的，只是要想象一张平面切不开一个球壳。

二、核技巧：把内积换掉

核 PCA 的全部机关只有一句话：PCA 只需要内积，那就把内积换成一个核函数。 PCA 的每一步都只用得到 X'X 这种内积形式。如果先把每个点映射成 φ(x)，那需要的是 φ(x)'φ(y)。 定义一个核函数 k(x, y) 直接给出这个内积，就再也不需要知道 φ 长什么样。 这就是「核技巧」：把「映射」这件事藏进内积里。 本节用的是高斯核 k(x, y) 等于 exp(-||x-y||²/2h²)。 它对应的是无穷维的特征空间——听起来很吓人，但代码里只多了一个距离矩阵。 写出来就是一串 exp，一行都不用循环。 流程四步：① 算 n×n 的核矩阵 K，K_ij 等于 k(x_i, x_j)；② 中心化 K； ③ 对中心化的 K 做特征分解；④ 第 i 个坐标等于第 i 个特征向量乘以对应特征值的平方根。 为什么第④步要乘特征值的平方根？因为核矩阵的特征值等于「特征空间里的方差」， 而我们要的是坐标本身，不是方差。开个根号就把尺度对上了。 这个细节很容易漏，漏了之后坐标会被压扁或者放大，虽然排序不变但尺度全错。

三、核矩阵必须中心化

这一节最容易踩的坑，也是代码上只差一行的坑：核矩阵必须中心化。 中心化的公式是 K 减去 1K/n 减去 K1/n 加上 1K1/n²，其中 1 是全 1 矩阵。 这四项目测很丑，但它的含义很朴素：在特征空间里减掉均值，跟普通 PCA 减均值是同一件事。 不中心化的后果有多严重？实测：第一主成分的最小值是 +0.0590、最大值 +0.7126，全部同号， 一个点都没落到负半轴。阈值取 0 相当于把 400 个点全判成一类，正确率正好 50.00%。 原因是：不中心化等于把特征空间的原点当成了数据中心。 同心圆整体离那个原点很远，每个点投影出来都是一大截正数，上面那点起伏全被这截常数淹掉。 这跟普通 PCA 不中心化的毛病是同一个，只是在核空间里更隐蔽——因为你根本看不见那个空间。 还有一个数值后果：不中心化时最大特征值是 94.0436，中心化后是 54.9875。 特征值变大不是「信息更多」，而是「多了一坨常数偏移」——这一坨会占据最大的那个特征方向， 把真正的结构挤到后面去。 工程上的检查办法很土但很有效：中心化后的核矩阵，特征值之和应该等于迹， 而且特征值里出现负数是正常的（核矩阵半正定但数值误差会让小的变负），取 max(0, λ) 就行。 但如果第一主成分全部同号，几乎可以断定是忘了中心化。

四、带宽是尺度，不是旋钮

高斯核只有一个参数 h，很多人把它当「超参数」，随便试几个取最好的。这个做法在观念上就错了。 h 决定的是「多远算远」：h 小，核值衰减快，每个点只认自己身边一小圈； h 大，衰减慢，所有点彼此都算「有点像」。 它必须和数据的实际尺度对上。

五、拿到的究竟是什么坐标

核 PCA 的第一主成分与半径的相关系数是 0.8128； 而线性 PCA 的第一主成分与半径的相关只有 0.0820。 这个对比说明：核 PCA 没做别的事，它只是把「半径」这个非线性坐标找出来了。 这也给了核方法一个很好用的直觉：核 PCA 找的是「在特征空间里方差最大」的方向，退回原始空间， 它对应的是某种「等值线弯弯的坐标」。 对同心圆，那个坐标就是半径；对别的数据，可能是角度、可能是距离某个中心的远近。 但要提醒一句：核 PCA 的两个坐标之间不一定「一个对应半径、一个对应角度」这种整齐的划分。 它只是在特征空间里按方差排队，退回来是什么样，取决于数据。 所以「第一主成分代表什么」这个问题在核方法里比在线性方法里更难回答。 还有一个绕不开的区别：核 PCA 没有显式的投影矩阵。 线性 PCA 里有一个 a，新点乘一下就有坐标； 核 PCA 里只有一个 n×n 的核矩阵，新点要和所有训练点算一遍核值才能得到坐标。 也就是说，新样本的坐标是「相对于训练集」算出来的，训练集变了坐标就变了。 工程上这一条很重要：如果模型要在线更新、或者要部署到资源受限的设备上， 得把训练集（或者它的一个子集）一起带上。

六、代价：n×n 的矩阵

核方法最贵的地方在于：不管原始维数是 2 还是 200，核矩阵永远是 n×n。 400 个点就是 400×400，一次全特征分解实测 36 到 39 毫秒。 复杂度大概是 n 的三次方。 n 从 400 涨到 4000，矩阵元素多 100 倍，特征分解要慢一千倍左右——几十毫秒会变成几十秒。 所以核方法在几万个点以上就不能硬算了。 常用出路有三条：一是 Nyström 近似，只算一部分列、剩下的用低秩近似补； 二是随机特征，用随机投影把核近似成有限维内积，退回线性 PCA； 三是先聚类或者抽样，把 n 压下来再算。 顺带对比一下：线性 PCA 的代价是 O(d³ 加 nd²)，d 是原始维数，和样本数 n 只是线性关系。 所以「数据量大、维数不高」用线性 PCA；「数据量中等、结构高度非线性」才考虑核 PCA。 选方法不是看哪个先进，是看哪个的瓶颈不在你身上。

七、动手 Demo：把同心圆撕开

这段代码把同心圆降到一维，先看线性 PCA 扫遍所有阈值能做到多少，再扫一遍核 PCA 的带宽， 最后验证「不中心化」这个坑有多致命。 注意第三张表：不中心化时第一主成分的最小值是正数——看到这个现象就该立刻想到「忘了减均值」。

import numpy as np

def dmat(A, B):
    return np.sqrt(np.maximum(((A[:, None, :] - B[None, :, :]) ** 2).sum(-1), 0.0))

def circles(n=200, seed=2):
    rg = np.random.default_rng(seed)
    a1 = rg.uniform(0.0, 2 * np.pi, n)
    a2 = rg.uniform(0.0, 2 * np.pi, n)
    X = np.vstack([np.column_stack([np.cos(a1), np.sin(a1)]),
                   np.column_stack([3 * np.cos(a2), 3 * np.sin(a2)])])
    y = np.r_[np.zeros(n), np.ones(n)]
    return X, y

def pca(X, k=1):
    Xc = X - X.mean(axis=0)
    C = Xc.T @ Xc / len(X)
    w, V = np.linalg.eigh(C)
    o = np.argsort(w)[::-1]
    return Xc @ V[:, o][:, :k]

def kpca(X, k=1, h=1.0, center=True):
    n = len(X)
    K = np.exp(-dmat(X, X) ** 2 / (2.0 * h * h))
    if center:
        one = np.full((n, n), 1.0 / n)
        K = K - one @ K - K @ one + one @ K @ one
    val, vec = np.linalg.eigh((K + K.T) / 2.0)
    o = np.argsort(val)[::-1]
    val, vec = val[o], vec[:, o]
    return vec[:, :k] * np.sqrt(np.maximum(val[:k], 0.0)), val

def thresh_acc(z, y):
    zz = np.sort(z.ravel())
    best = 0.0
    for c in np.unique((zz[1:] + zz[:-1]) / 2.0):
        a = float(np.mean((z.ravel() > c).astype(float) == y))
        best = max(best, a, 1.0 - a)
    return best

X, y = circles(200, 2)
Z1 = pca(X, 1)
print('线性PCA 最优阈值正确率 %.4f' % thresh_acc(Z1, y))

for h in (0.5, 1.0, 1.5, 2.0, 3.0):
    Zk, v = kpca(X, 1, h)
    a0 = float(np.mean((Zk[:, 0] > 0).astype(float) == y))
    print('h=%.1f 阈值0正确率 %.4f 最大特征值 %.4f'
          % (h, max(a0, 1 - a0), v[0]))

Zc, vc = kpca(X, 1, 1.0)
Zu, vu = kpca(X, 1, 1.0, center=False)
print('不中心化 最小 %+.4f 最大 %+.4f 最大特征值 %.4f'
      % (Zu[:, 0].min(), Zu[:, 0].max(), vu[0]))
print('中心化   最小 %+.4f 最大 %+.4f 最大特征值 %.4f'
      % (Zc[:, 0].min(), Zc[:, 0].max(), vc[0]))

r = np.sqrt((X ** 2).sum(1))
print('线性PC1与半径相关 %+.4f' % np.corrcoef(r, Z1[:, 0])[0, 1])
print('核PC1与半径相关   %+.4f' % np.corrcoef(r, kpca(X, 2, 1.5)[0][:, 0])[0, 1])

实跑输出：

--------------------------------------------------------------------
同心圆降到一维：PCA 切到 72.5%，核 PCA 换对带宽就到 100%
--------------------------------------------------------------------
  数据                     400 个二维点：内圈半径 1 共 200 点，外圈半径 3 共 200 点
  任务                     只用一维坐标把两圈分开

  先看线性投影能做到什么程度。PCA 降到一维，然后把所有可能的阈值都试一遍：
  最好的一刀               阈值 0.9078 时正确率 0.7250
  一维上的 1 近邻          0.8500
  **一条直线穿过圆心，内圈被压到中间、外圈被压到两头，中间那段必然叠在一起。**
  任何一刀都切不干净，这不是阈值没调好，是投影本身把信息丢了。

  核 PCA 换个思路：先用高斯核 exp(-d²/2h²) 把点映射到无穷维，
  在那个空间里做 PCA，再把结果拉回一维。唯一要调的参数是带宽 h。
  带宽 h   阈值 0 的正确率   最大特征值   内圈均值   外圈均值
  ----------------------------------------------------------
     0.5            0.7650      37.9424    -0.0079    +0.0079
     1.0            1.0000      54.9875    -0.3618    +0.3618
     1.5            0.9775      59.7835    -0.3142    +0.3142
     2.0            0.5775      61.6964    +0.1049    -0.1049
     3.0            0.5450      55.1299    +0.0469    -0.0469
  **h 等于 1.0 的时候两圈被彻底撕开：内圈全在负半轴，外圈全在正半轴，阈值取 0 就是完美分界。**
  h 太小（0.5）时核值只在自己身上还算大，点和点之间连不起来，图成了一张碎石子；
  h 太大（2.0 以上）所有点和所有点都差不多，核矩阵拉平，又退回到线性投影的水平。
  **带宽不是那种「调调看」的超参数，它必须和数据本身的尺度对上——这里内圈半径就是 1，h 取 1.0 刚好。**

  再验一件事：第一主成分是不是真的抓住了半径？
  第一主成分来自    与半径的相关
  -----------------------------
  线性 PCA               +0.0820
  核 PCA（h=1.5）        +0.8128
  **线性 PCA 的第一主成分只是一条穿过圆心的直线方向，和半径几乎不相关；核 PCA 的第一主成分本身就是半径的单调函数。**
  也就是说，核 PCA 没做别的事，它只是把「半径」这个非线性坐标找出来了。

  最后一个坑，也是新手最常踩的：核矩阵必须中心化。
  核矩阵     最大特征值   第一主成分的最小值    最大值   阈值 0 的正确率
  ---------------------------------------------------------------------
  不中心化      94.0436              +0.0590   +0.7126            0.5000
  中心化        54.9875              -0.4570   +0.5069            1.0000
  **看后两列：不中心化时第一主成分全部同号，一个点都没落到负半轴，阈值 0 相当于把所有人都判成一类。**
  **原因是不中心化等于把特征空间的原点当成了数据中心。同心圆整体离那个原点很远，于是每个点投影出来都是一大截正数，上面那点起伏全被淹掉。**
  **中心化做的事就是「在特征空间里减掉均值」——核矩阵那一串 K - 1K/n - K1/n + 1K1/n² 全是这一步展开的结果。**
  这个错误在代码上只差一行，在现象上表现为「结果看起来完全是随机的」。

  代价也顺手量一下（跑 20 次取最快的一次，避开机器忙碌的干扰）：
  n=400 一次核 PCA         37.8 毫秒
  核矩阵是 n×n 的，要做一次全特征分解，复杂度大约 n 的三次方。
  **n 从 400 涨到 4000，矩阵元素多 100 倍，特征分解要慢一千倍左右。**
  所以核方法在几万个点以上就不能硬算了，得用 Nyström 之类的近似。

第一行就是线性方法的判决：扫遍所有阈值，最好 0.7250。 注意这是「扫遍」，不是随手取 0——已经把阈值调到最优了还是切不干净，说明问题出在投影本身， 不在阈值。 带宽那张表里，h 等于 1.0 那一行是 1.0000，阈值取 0 就是完美分界； 而 h 取 0.5 时内圈均值和外圈均值只差 0.0158，几乎重合。 这一行最能说明「带宽是尺度」：内圈半径是 1，h 取 1.0 的时候核的衰减和几何尺度刚好对上。 h 从 1.0 到 2.0，正确率从 1.0000 掉到 0.5775，中间没有平缓的过渡。 所以实际用的时候，h 要在对数尺度上扫（0.01、0.03、0.1、0.3、1、3……），不能线性地试几个整数。 后面两组数里，最大的信息量在不中心化那一行的符号上：最小值 +0.0590、最大值 +0.7126， 全是正的。 第一主成分全部同号，就意味着「阈值 0」这个判据失效——它会把所有人判成一类，正确率 0.5000。 这不是核 PCA 的错，是漏了一步中心化。看到「第一主成分全同号」就该往这里查。 最后两行把核 PCA 干的事说白了：线性 PC1 与半径的相关只有 0.0820， 核 PC1 与半径的相关是 0.8128。 核 PCA 不是找到了什么神秘的东西，它就是把「半径」这个弯弯的坐标翻出来了。

30 秒实验

1. 把 circles 的两个半径从 1 和 3 改成 1 和 1.2（两圈挨得很近）。 核 PCA 的最优 h 会明显变小——还是那句话，h 要跟几何尺度对上。
2. 把 kpca 里的 np.sqrt(np.maximum(val[:k], 0.0)) 改成 val[:k]（不乘平方根）。 坐标会被整体放大或压扁，但排序不变、阈值 0 的正确率也不变。想清楚： 什么时候尺度会真的影响结果？（提示：算近邻、算距离的时候。）
3. 把数据集换成 swiss(600, 1) 那种三维数据，用核 PCA 降到二维，再算一次 10 近邻保持率。 会发现核 PCA 在上面并不比线性 PCA 好多少——核 PCA 也不是万能的，它仍然是一种「全局」方法。

常见误区

1. 「核 PCA 就是 PCA 加了核，肯定更好」。它换的是特征空间，代价是多了带宽这个必须调对的参数， 以及 n×n 的内存与时间。在这份同心圆上它 100% 对线性 PCA 的 72.50%； 但在上一节的瑞士卷上，它的近邻保持率并没有明显超过线性 PCA。
2. 「忘了中心化只是差一点点」。差的是全部。 实测不中心化时第一主成分全部同号，正确率从 1.0000 掉到 0.5000，等于完全失效。 最省事的检查就是看第一主成分是不是全同号。
3. 「h 越小越精细，效果越好」。 h 太小，核矩阵接近对角阵，每个点只认得自己， 什么结构都学不到——实测 h 取 0.5 时两圈的均值差只有 0.0158。 h 太大，核矩阵被拉平，退化成线性方法。两头都不行。
4. 「核 PCA 得到的坐标可以直接和原始坐标对照解释」。 它是在无穷维特征空间里算出来的，坐标的绝对值和符号都没有物理含义， 只有「相对的远近」有意义。 说要解释，只能说明「第一主成分与某个原始量的相关是多少」这类间接关系。

核 PCA 走的是「先升维、再线性」这条迂回路。 下一节走另一条：不升维也不投影，只利用那张邻接图本身——把数据当成一张弹簧网，让它自己松开。

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

相关笔记

• 主成分分析
• 核模型
• 拉普拉斯特征映射
• 🤖 图解机器学习 · 索引