局部线性嵌入

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

上一章末尾讲 LLE 时说过：核心思路好懂，但每个公式的具体推导「需要下节课再跟进」——这节课就来兑现这句话。我们把 LLE（Local Linear Embedding，局部线性嵌入）的两步数学完整推一遍：第一步在高维空间里学「每个点由邻居线性重构」的权重 $W$，第二步带着这份 $W$ 去低维空间里找一组坐标 $Y$——最后你会发现，兜兜转转，落脚点又是我们的老朋友：特征值与特征向量。

💡 核心导读

• 嵌入（embedding）：把不容易直接表示的量（词、图、情感、句子……）变成空间里的一个向量——机器学习里特别常用的说法；LLE 是它在降维领域的一个具体实现。
• LLE 的核心假设：高维空间里每个点，都可以由它最近的 $k$ 个邻居线性组合出来——「我是谁，由我身边的朋友定义」。
• 两步走：第一步对每个点求重构权重 $W$（带约束的优化，拉格朗日乘子出场）；第二步固定 $W$、求低维坐标 $Y$（同样带约束，同样拉格朗日）。
• 落脚点是特征分解：第二步化成 $My=\mu y$ 的形式——$M=(I-W)^T(I-W)$，取最小的非零特征值对应的特征向量做嵌入（丢弃的那个零特征值，恰好就是平移方向）。
• 一个容易混的点：LLE 名字里的 Linear 指局部线性重构，不是说整体映射是「乘一个矩阵」的线性降维——$Y$ 不是拿 $X$ 乘出来的。

嵌入：把不容易表示的东西变成向量

先说名字。embedding 中文一般叫嵌入。大家学到深度学习的时候会碰到 Google 的 word2vec：把任何一个词映射成一个向量，中文叫词嵌入。后来你会发现有各种各样的嵌入办法——所谓的嵌入，就是把一些不太容易表示的量，变成空间里的一个向量。图可以嵌入、词可以嵌入、情感可以嵌入、句子也可以嵌入——讲师的原话是「embedding everything, everything can be embedded」。在机器学习里边，「嵌入」这个说法特别常用。

我们这里讲的是 Local Linear Embedding：通过局部线性重构的方式，在低维空间里重新安放数据，来满足它们在高维空间里的特性。

LLE 的直观：我是谁，由我的邻居定义

在高维空间 $\mathbb{R}^N$ 里有很多数据。挑出我们需要研究的一个点 $x_i$（讲师板书上的红点）：我们认为这个红点可以表示成它周围几个最近邻的线性组合。周围那几个点（板书上的绿点）是 $\eta_1, \eta_2, \ldots, \eta_k$，每个绿点连向红点的那条线，就带一个权重 $w$：

$$x_i \approx \sum_{j \in \text{邻居}} w_j \, \eta_j$$

我们希望「$x_i$ 减去它的重构」的平方最小，而权重满足一个条件：周围这 $k$ 个点的权重加和等于 1。

讲师给了一个好玩的类比。我是一个陌生人，我想知道我的数学成绩是多少——很惊奇的是，我可以去找我最好的几个朋友：我经常跟谁在一起玩，就把他们的数学成绩做一个加权平均，那就是我的数学成绩。因为我们相信「人以类聚，物以群分」：学习好的跟学习好的小朋友在一起，学习差的跟学习差的在一起（当然这不绝对，就是个好玩的例子）。这样即使我并不知道你的数学成绩，也可以通过你经常一起玩的小朋友的数学成绩线性组合出来。这个线性组合的 $w$ 怎么组合？——这正是我们要学习的部分：它就是我们的参数、我们的权重。

再往降维的方向问一句。刚才研究的也许不是数学成绩一门，而是数理化三门成绩——那是一个三维的表示；映射到一个低维空间，可能就变成了理科平均成绩（一维）。那我们想知道：一个人的理科平均成绩，和他周围最好小伙伴的理科平均成绩之间，还是同样的线性组合关系吗？如果希望这个关系在低维空间里也保持，低维的表示 $Y$ 应该怎么取？——「直观上我们会问：是不是像 PCA 那样，$X$ 乘上一个矩阵 $A$ 就映射到 $Y$ 空间了？」推导到最后你会发现，答案比「乘一个矩阵」更精巧：$Y$ 直接取自一个矩阵的特征向量——这也正是 LLE 属于非线性方法的地方。下面就把这套数学完整推出来。

第一步：求重构权重

从 $e_i$ 到矩阵形式

先定义每个点的误差 $e_i$。给定 $x_i$（$i$ 从 1 到 $N$，共 $N$ 个数据），它需要对周围的几个点构造出来——周围那几个最近邻是 $\eta_j$，对应的权重是 $w_j$。我们要最小化的目标（对 $j$ 从 1 到 $k$ 加和）：

$$e_i = \Big\Vert x_i - \sum_{j=1}^{k} w_j \, \eta_j \Big\Vert^2, \qquad \sum_{j=1}^{k} w_j = 1$$

推导的第一个小技巧，是把 $x_i$ 自己「塞进」加和里。看这个变形（讲师现场一步步带大家验证过相等）：

$$x_i - \sum_j w_j \eta_j = \sum_j w_j x_i - \sum_j w_j \eta_j = \sum_j w_j \left( x_i - \eta_j \right)$$

第一步成立是因为 $\sum_j w_j = 1$，所以 $\sum_j w_j x_i = x_i$。于是误差改写为 $e_i = \big\Vert \sum_j w_j (x_i - \eta_j) \big\Vert^2$——把 $w_j$ 从中间「挖」了出来，剩下的是 $x_i$ 与每个邻居的差值向量。

把 $k$ 个差值向量按列排成矩阵 $X$（每列是 $x_i - \eta_j$，$j = 1, \ldots, k$），权重排成向量 $w = (w_1, \ldots, w_k)^T$，则 $e_i = \Vert Xw \Vert^2 = w^T (X^T X) w$。注意 $X$ 里边的 $\eta_j$ 本身也是 $x$（周围的数据点），所以这一步只是换个写法。目标变成：

$$\min_w \; w^T (X^T X) w \qquad \text{s.t.} \qquad \mathbf{1}^T w = 1$$

这里 $\mathbf{1}$ 是全 1 列向量（$(1,1,\ldots,1)^T$）——$\mathbf{1}^T w$ 就是 $w_1 + w_2 + \cdots + w_k$。（讲师口播时把它说成「单位阵」，按写法是全 1 向量；「每一行加和为 1」和「$\mathbf{1}^T w = 1$」是同一件事的两种说法。）

拉格朗日乘子：老朋友又出场了

这是一个有约束的优化问题——大家应该已经很熟悉了：用拉格朗日乘子。构造拉格朗日函数（把约束改写成 $1 - \mathbf{1}^T w = 0$）：

$$L(w, \lambda) = w^T (X^T X) w + \lambda \left( 1 - \mathbf{1}^T w \right)$$

讲师在这里专门回扣了一句：「如果忘了，可以去看 SVM 前面那部分的 note 和课程——我花了几乎一节课去讲拉格朗日乘子、讲有条件优化」（对应本课第 10 章的拉格朗日乘子法三节：无约束到等式约束、平行条件、KKT）。

对 $w$ 求导并令其为零。这里讲师停下来讲了一个很多人不清楚的细节——向量求导：一个数对一个向量求导，相当于这个数对向量的每一个分量分别求导，所以结果还是一个向量。$w$ 是 $k$ 维的（$w_1$ 到 $w_k$，对应周围最近的 $k$ 个点），于是：

$$2 (X^T X) w = \lambda \mathbf{1} \quad \Rightarrow \quad w = \frac{\lambda}{2} (X^T X)^{-1} \mathbf{1}$$

对 $\lambda$ 求导则原样回收约束 $1 - \mathbf{1}^T w = 0$，没有新信息。

归一化收尾：$\lambda$ 自己约掉了

代入约束 $\mathbf{1}^T w = 1$ 可解出 $\lambda/2$；但更省事的做法是直接「先解再归一」：$w$ 反正正比于 $(X^T X)^{-1} \mathbf{1}$，比例常数是多少不影响「加和为 1」——那就把每个分量算出来之后，底下做个 normalization（除以总和）就完了。讲师给中间那个矩阵起了个名字：$Q = X^T X$（$k \times k$；它正是那 $k$ 个差值向量的局部 Gram 矩阵，$X^T X / n$ 的口径就是协方差矩阵）。于是：

$$w_j = \frac{\left( Q^{-1} \mathbf{1} \right)_j}{\mathbf{1}^T Q^{-1} \mathbf{1}}, \qquad j = 1, \ldots, k$$

（讲师现场命名 $Q$ 时绕了一下——先说成「$X^T X$ 的负一次幂叫 $Q$」、后又写 $Q^{-1}$，按最终用法，$Q$ 就是 $X^T X$ 本身，权重用 $Q^{-1}$。）注意 $x$ 固定了，找最近邻是确定的，$\eta_j$ 固定，$Q$ 就完全确定——通过 $x$ 的关系，就可以找到 $w$ 的关系，整个链条是干净的线性代数。

用代码把这一步真正算一遍。讲义设定：平面上三个点，用两个邻居 $\eta_1=(0,0)$、$\eta_2=(4,0)$ 重构 $x_0=(1, 0.5)$（刻意让 $x_0$ 不在两邻居的连线上）：

import math

eta1 = (0.0, 0.0)
eta2 = (4.0, 0.0)
x0 = (1.0, 0.5)  # 讲义设定: 待重构的点, 不在两邻居连线上

Z = [[x0[d] - eta1[d] for d in range(2)],
     [x0[d] - eta2[d] for d in range(2)]]          # 每行一个差值向量 z_j
Q = [[sum(Z[r][d] * Z[c][d] for d in range(2)) for c in range(2)] for r in range(2)]

det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0]
u = [(Q[1][1] - Q[0][1]) / det, (Q[0][0] - Q[1][0]) / det]   # u = Q^{-1} 1
w = [ui / sum(u) for ui in u]

print("差值向量 z1, z2 =", Z)
print("局部 Gram 矩阵 Q =", [[round(v, 4) for v in row] for row in Q])
print("归一化前 u = Q^{-1} 1 =", [round(ui, 4) for ui in u])
print("重构权重 w =", [round(wi, 4) for wi in w], " 和 =", round(sum(w), 4))
recon = [w[0] * eta1[d] + w[1] * eta2[d] for d in range(2)]
print("重构点 =", [round(v, 4) for v in recon], " 重构误差 =", round(math.dist(recon, x0), 4))

差值向量 z1, z2 = [[1.0, 0.5], [-3.0, 0.5]]
局部 Gram 矩阵 Q = [[1.25, -2.75], [-2.75, 9.25]]
归一化前 u = Q^{-1} 1 = [3.0, 1.0]
重构权重 w = [0.75, 0.25]  和 = 1.0
重构点 = [1.0, 0.0]  重构误差 = 0.5

三个看点：权重是 $(0.75, 0.25)$——离谁近，谁的权重就大（$x_0$ 的投影离 $\eta_1$ 近，$\eta_1$ 分到 0.75）；权重加和恰好 1.0；重构点落在 $(1.0, 0.0)$——正是 $x_0$ 在「邻居连线」上的投影，剩余误差 0.5 是它离开那张「局部线性片」的垂直距离。约束下的最优重构，不一定是零误差——最优的权重给出的是邻居所能拼出的、离 $x_0$ 最近的那个点。

第二步：求低维嵌入

把整个问题搬到 $Y$ 空间

有了每个点的权重，把高维数据映射到低维：$\mathbb{R}^N \to \mathbb{R}^M$（$M$ 远小于 $N$，这才起到降维的作用）。构造新空间里的误差——形式跟刚才一模一样，只是 $x$ 换成了 $y$：

$$E_Y = \sum_{i=1}^{N} \Big\Vert y_i - \sum_{j} w_{ij} \, y_j \Big\Vert^2$$

$w_{ij}$ 就是刚才在高维空间里辛苦求出来的那份权重（$j$ 不是 $i$ 的邻居时为 0）；我们要找一组 $y$，让同样的 $w$ 在低维空间里照样把每个点拼出来。

推导用跟第一步平行的「挖中间项」技巧：把 $y$ 从减法中间提出来，$y_i - \sum_j w_{ij} y_j = \sum_j (\delta_{ij} - w_{ij}) y_j$，平方展开成转置乘积（讲师边写边说「我们刚才所做的事情是一样的」，中途还停下来重誊了一次板书）。于是整个目标写成矩阵形式，中间那坨东西单独命名——这就是 $M$ 矩阵：

$$M = (I - W)^T (I - W), \qquad E_Y \propto \; \text{tr}(Y^T M Y)$$

两个限制条件：位置可调，所以钉死

到这步还没完——$Y$ 还得加两条限制，否则解不唯一。讲师给的动机非常直观：

第一条：中心化，$\sum_i y_i = 0$。 把高维数据映射到低维空间时，它的位置是可调的——映射到以 0 为中心、以 1 为中心、以 100 为中心都可以，全体的平移不影响「邻居关系」这个我们真正关心的东西。既然可调，那就干脆强制：所有 $y_i$ 加和等于零，以零为中心排列。

第二条：单位协方差，$Y^T Y / n = I$。 让低维表示的协方差矩阵是单位阵——各维度去相关、统一尺度，不偏袒任何一个低维方向。

又是拉格朗日，又是特征分解

优化 $E_Y$、给定这两条限制，又可以构造拉格朗日函数了（第二个了——这门课走到这里，拉格朗日乘子几乎每一章都要出场一次）。对 $y$ 求导令其为零，得到的是一个特别眼熟的形式：

$$M y = \mu y$$

（$\mu$ 是新的拉格朗日系数，讲师随口起名叫「缪」；左边系数、约分合并的过程跟第一步完全同款。）——这不就是特征分解吗！乘上 $M$ 这个矩阵，相当于乘上一个常数：我们想求的，就是 $M$ 矩阵的特征值和特征向量。

所以对 $M$ 做特征分解，然后挑特征向量。这里有一个方向性的坑，值得单独点出来：讲师口播时三次说成「找最大的几个特征值、保留最大的几个特征向量」——这是说反了。我们全程都在最小化误差 $E_Y = \text{tr}(Y^T M Y)$（Rayleigh 商的最小化方向），标准解法是取最小的特征值；而 $M = (I-W)^T(I-W)$ 又恰好有一个天然的零特征值——每一行的权重加和为 1，意味着 $(I-W)\mathbf{1} = 0$，也就是 $M \mathbf{1} = 0$：全 1 向量是零特征值的特征向量，它就是被第一条约束排除掉的平移方向。所以最终取法是：丢弃最小的零特征值，取第 2 到第 $M+1$ 小的特征值所对应的特征向量，按行排成低维嵌入 $Y$。

用代码把这两件事一起看到：数据取半圆弧上的 9 个点（一条嵌在二维里的一维流形，正好回扣上一节「流形」），跑完整 LLE 两步——

import math

n = 9
arc = [(math.cos(math.pi * i / (n - 1)), math.sin(math.pi * i / (n - 1))) for i in range(n)]

def dist(a, b):
    return math.hypot(a[0] - b[0], a[1] - b[1])

# 第一步: 每个点取 2 个最近邻, 求 w ∝ Q^{-1} 1 再归一化
W = [[0.0] * n for _ in range(n)]
for i in range(n):
    order = sorted((j for j in range(n) if j != i), key=lambda j: dist(arc[i], arc[j]))
    nbrs = order[:2]
    Z = [[arc[i][d] - arc[j][d] for d in range(2)] for j in nbrs]
    Q = [[sum(Z[r][d] * Z[c][d] for d in range(2)) for c in range(2)] for r in range(2)]
    det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0]
    if abs(det) < 1e-12:  # 近邻几乎共线时兜底均分
        w = [0.5, 0.5]
    else:
        u = [(Q[1][1] - Q[0][1]) / det, (Q[0][0] - Q[1][0]) / det]
        w = [ui / sum(u) for ui in u]
    for idx, j in enumerate(nbrs):
        W[i][j] = w[idx]

# 第二步: M = (I - W)^T (I - W)
M = [[sum(((1.0 if a == b else 0.0) - W[a][b]) * ((1.0 if a == c else 0.0) - W[a][c]) for a in range(n)) for c in range(n)] for b in range(n)]

def jacobi_eig(A, tol=1e-12, max_sweeps=100):
    """雅可比迭代: 对称矩阵特征分解, 返回(特征值, 特征向量按列)"""
    n = len(A)
    a = [row[:] for row in A]
    V = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
    for _ in range(max_sweeps):
        off = math.sqrt(sum(a[p][q]**2 for p in range(n) for q in range(n) if p != q))
        if off < tol:
            break
        for p in range(n - 1):
            for q in range(p + 1, n):
                if abs(a[p][q]) < 1e-15:
                    continue
                theta = 0.5 * math.atan2(2 * a[p][q], a[q][q] - a[p][p])
                c, s = math.cos(theta), math.sin(theta)
                for k in range(n):
                    akp, akq = a[k][p], a[k][q]
                    a[k][p] = c * akp - s * akq
                    a[k][q] = s * akp + c * akq
                for k in range(n):
                    apk, aqk = a[p][k], a[q][k]
                    a[p][k] = c * apk - s * aqk
                    a[q][k] = s * apk + c * aqk
                for k in range(n):
                    vkp, vkq = V[k][p], V[k][q]
                    V[k][p] = c * vkp - s * vkq
                    V[k][q] = s * vkp + c * vkq
    return [a[i][i] for i in range(n)], V

vals, V = jacobi_eig(M)
order = sorted(range(n), key=lambda i: vals[i])
print("M 的特征值(升序):", [f"{vals[i]:.6f}" for i in order])
v0 = [V[i][order[0]] for i in range(n)]
print("最小特征值对应向量(归一化后, 应为常数=平移方向):", [round(v / max(abs(x) for x in v0), 3) for v in v0])
y = [V[i][order[1]] for i in range(n)]
print("1D 嵌入 y(最小非零特征值对应的特征向量):", [round(v, 3) for v in y])
print("弧上顺序 0..8 对应 y 是否单调:", all(y[i] < y[i+1] for i in range(n-1)))

M 的特征值(升序): ['0.000000', '0.000039', '0.020896', '0.167052', '0.671359', '1.755549', '3.206505', '7.882226', '7.906106']
最小特征值对应向量(归一化后, 应为常数=平移方向): [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
1D 嵌入 y(最小非零特征值对应的特征向量): [-0.512, -0.39, -0.261, -0.13, 0.0, 0.13, 0.261, 0.39, 0.512]
弧上顺序 0..8 对应 y 是否单调: True

三行输出对应三件事：最小的特征值是零（数值上 0.000000），它对应的特征向量归一化后是清一色的 1.0——正是全 1 平移方向，扔掉；最小非零特征值（0.000039，一维流形的展开模式在这个谱里能量极低）对应的特征向量，就是 9 个点的 1D 嵌入——从 -0.512 到 0.512 近似等距地单调排开：半圆弧被「拉直」成了一条线段，而且点距保持得很好。一条弯在二维里的一维流形，就这么被 LLE 展开了——这就是「保持局部线性关系」换来的全局效果。

全流程总结：两个问题，一个特征分解

最后把整件事串一遍。LLE 是一个从高维空间映射到低维空间的运算，满足的条件是：在高维空间里，每一个点由它周围的邻居线性组合构造；到了低维空间，这批点（变成 $y$ 之后）仍然通过同一份 $w$ 组合到一起。所以问题归结为两个：

1. $W$ 怎么估计？——刚才推的：每个点找 $k$ 个最近邻，解带 $\mathbf{1}^T w = 1$ 约束的最小化，$w \propto Q^{-1}\mathbf{1}$ 再归一化（$Q = X^TX$ 完全由数据决定）。
2. $X$ 怎么映射到 $Y$？——构造 $M = (I-W)^T(I-W)$，对 $M$ 做特征分解，取最小非零的几个特征值对应的特征向量。比如 $M$ 一共是十维的，有十个特征值；保留（最小的非零的）三个特征向量，就相当于把 10 维映射到了 3 维。

还有一个容易混的点再钉一次：讲师总结时说「$X$ 映射到 $y$ 是一个降维，是个线性的降维，一定是乘一个矩阵」——LLE 名字里的 Linear，指的是局部线性重构，不是全局线性映射。$Y$ 的每一行是 $M$ 的特征向量按点排布，不存在 $X \cdot A = Y$ 这样的变换矩阵（那是 PCA 的形态）；LLE 直接输出每个点的低维坐标，这正是它「非线性」的所在。

📝 动手练一练

1. 概念辨析：PCA 和 LLE 都把落脚点放在「特征值与特征向量」上，但方向完全相反。请从三个方面对比两者：各自假设数据有什么样的结构？各自保留什么？解的形态差在哪（有没有 $X \cdot A = Y$ 形式的映射）？

2. 代码实验：把上面块1 代码里的 $x_0$ 从 $(1, 0.5)$ 改成 $(5, 0.5)$——跑到 $\eta_2$ 的外侧去了。先猜权重会变成什么样，再跑代码验证。这说明 LLE 的权重有什么性质？

# 为保证本块可独立运行, 重写权重求解
import math

eta1 = (0.0, 0.0)
eta2 = (4.0, 0.0)
x0 = (5.0, 0.5)  # 跑到 eta2 外侧

Z = [[x0[d] - eta1[d] for d in range(2)],
     [x0[d] - eta2[d] for d in range(2)]]
Q = [[sum(Z[r][d] * Z[c][d] for d in range(2)) for c in range(2)] for r in range(2)]

det = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0]
u = [(Q[1][1] - Q[0][1]) / det, (Q[0][0] - Q[1][0]) / det]
w = [ui / sum(u) for ui in u]

print("差值向量 z1, z2 =", Z)
print("局部 Gram 矩阵 Q =", [[round(v, 4) for v in row] for row in Q])
print("归一化前 u = Q^{-1} 1 =", [round(ui, 4) for ui in u])
print("重构权重 w =", [round(wi, 4) for wi in w], " 和 =", round(sum(w), 4))
recon = [w[0] * eta1[d] + w[1] * eta2[d] for d in range(2)]
print("重构点 =", [round(v, 4) for v in recon], " 重构误差 =", round(math.dist(recon, x0), 4))

差值向量 z1, z2 = [[5.0, 0.5], [1.0, 0.5]]
局部 Gram 矩阵 Q = [[25.25, 5.25], [5.25, 1.25]]
归一化前 u = Q^{-1} 1 = [-1.0, 5.0]
重构权重 w = [-0.25, 1.25]  和 = 1.0
重构点 = [5.0, 0.0]  重构误差 = 0.5

本章小结

这一节把 LLE 的数学完整推完了。直观：每个点由它最近的 $k$ 个邻居线性重构——「我是谁，由我身边的朋友定义」；权重加和为 1，是参数、是学出来的。第一步：对每个点解带约束的最小化 $w^T (X^TX) w$——拉格朗日乘子（第 10 章的老朋友）给出 $w \propto Q^{-1}\mathbf{1}$，归一化收尾；$Q$ 完全由数据决定，链条干净。第二步：把同样的重构关系搬到低维，误差写成 $\text{tr}(Y^T M Y)$（$M = (I-W)^T(I-W)$）；两条限制（以零为中心、单位协方差）钉死平移和尺度自由度；再上一次拉格朗日，得到 $My = \mu y$——特征分解。取最小非零特征值对应的特征向量（全 1 的零特征向量是平移方向，丢弃），半圆弧的实验亲眼看到：弯曲的一维流形被拉直成近似等距的线段。最后记住那个容易混的点：LLE 的 Linear 是局部线性，$Y$ 不是 $X$ 乘矩阵乘出来的——这正是它与 PCA 的分野。

📋 行动清单

• 复述两步：合上讲义，说出 LLE 两步各自在优化什么、约束各是什么、各自用什么数学工具解。
• 验证 $M\mathbf{1}=0$：在块2 代码里加一行，直接计算 $M \cdot \mathbf{1}$，亲眼看到它等于零向量。
• 改近邻数：把块2 的近邻数从 2 改成 3（注意 $k=3$ 时三个邻居的 $Q$ 可能接近奇异——体会讲义代码里那个兜底分支为什么存在）。
• 对照 PCA：把 16.3 节 PCA 的「取最大特征值」与本节「取最小非零特征值」放在一起，用自己的话讲清为什么方向相反。

—— 小象教研组