SVD 降维

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

上一节的 PCA 用「乘一个矩阵」把高维压成低维。这一节讲另一种形式——它跟 PCA 类似，但效果在某种程度上并不输于它：奇异值分解（SVD）。我们还用它解一个自然语言里的经典难题：一个词有两副面孔（rock 是石头还是摇滚？），机器怎么分辨？最后，把上次课没讲完的非线性降维（Isomap 与 LLE）一并收线。

💡 核心导读

• SVD 的形式：任何 $n \times d$ 矩阵都能分解成三个矩阵的乘积 $U \Sigma V^T$；中间的 $\Sigma$ 是对角阵，对角线上的值叫奇异值。
• 降维机制：把 $\Sigma$ 里的值从大到小重排，取最大的 $k$ 个、其余当作 0——剩下三个「瘦身」矩阵的乘积就是原矩阵的逼近 $\hat{A}$。
• 潜在语义分析（LSA/LSI）：对词-文档矩阵做 SVD，文档被压进低维「概念空间」——语义的距离可以用余弦夹角来量。
• 三方法递进：LSA 纯代数 → PLSA 加概率形式和一个随机变量 → topic model（LDA）再假设一层随机变量——第 15 章三件套的代数侧起点。
• 非线性降维：乘矩阵永远是线性变换；面对瑞士卷、麻花状的流形数据，Isomap 把距离换成图上最短路（Dijkstra），LLE 用「邻居的线性组合」约束映射。

SVD：三个矩阵的乘积

奇异值分解（Singular Value Decomposition）做的事，是把一个完整的 $n \times d$ 矩阵，分解成三个矩阵的乘积：

$$\underbrace{A}_{n \times d} = \underbrace{U}_{n \times r} \cdot \underbrace{\Sigma}_{r \times r} \cdot \underbrace{V^T}_{r \times d}$$

三个矩阵乘起来，可以完整地把原矩阵表述出来。中间的 $\Sigma$ 是一个对角阵——只有对角线上有值，其他地方全为 0（讲师讲到这里的记号时先说成「欧米伽」，马上自我纠正：是西格玛 $\Sigma$）。这个对角阵非常关键：$\Sigma$ 对角线上的值，就相当于 PCA 里特征值的角色——我们要选取最大的几个，其余的去掉。SVD 做降维时干的就是同一件事：把 $\Sigma$ 里的值重新排列，使最大的值排在左上角（$\sigma_1, \sigma_2, \ldots$，如同 $\lambda_1, \lambda_2, \ldots$），然后取最大的几个值。

实现上大家完全不用担心：Python、MATLAB、R 任何一个语言里都有现成的——PCA 和 SVD 都是代数里传统的运算，你不用写自己的程序，打一个 svd 它就自动分解，打 pca 也可以分解。

那降维具体怎么发生？看这张「结构图」：最开的时候矩阵是 $V \times M$ 维（词数 $\times$ 文档数），分解成 $V \times R$、$R \times R$、$R \times M$ 三个矩阵的乘积，完全等价于 $A$，没有任何问题。如果想降到 $K$ 维，就在 $\Sigma$ 里选最大的 $k$ 维保留，把剩下的部分全部删掉——那些位置默认为 0，而灰色部分乘以 0 就全部消失了。最后剩下的三个「黑色长条」：$V \times k$、$k \times k$、$k \times M$，乘起来仍然是一个 $V \times M$ 的矩阵，记作 $\hat{A}$——它是原始矩阵 $A$ 的一个逼近，去掉的只是无关紧要的部分。

LSI：rock 是石头还是摇滚

现在看一个 classical（经典）的问题，看 SVD 怎么在实际中用。给定一个 term by document matrix（词-文档矩阵）。它背后的假设我们在第 15 章就说过：自然语言处理里，我们假设任何词之间不存在前后关系、不存在顺序上的关联——每一个文本就是一组词的袋子（a bag of words，词袋模型）。

矩阵里的数是词在文档里出现的频率。看讲师的例子（词典里六个词：rock、granite、marble、music、song、band）：

• 第一个文档：rock 出现了 2 次，granite 1 次，marble 1 次；music、song、band 没有出现。（讲师口播「granite 是石灰石」有口误——granite 是花岗岩，limestone 才是石灰石；marble 是大理石，rock 是石头。）
• 第四个文档：rock 也出现了 2 次，但没有 granite 和 marble——出现的是 music 和 song。

所以问题大家已经很清楚了：rock 这个词有两副非常清晰的面孔——一类是 rock music（摇滚乐），一类是真正的石头。这正是自然语言里非常重要的问题：怎么做 disambiguation（歧义消除）？因为 rock 本身有不同的 semantic meaning（语义信息）。我们希望的是：给了一个新文本，能不能看出它里面的 rock 在描述什么？这个新文本跟语料库里的哪篇文章距离最近？

讲师的例子里有 6 个训练文本，矩阵的每个值代表词出现的频率。现在给一个新的 query（查询）——「rock marble」，其他词没有出现——问这个 query 跟 6 个文档之间，哪一个距离最近？

以前研究机器学习的人都知道，语义是非常难处理的情况——很长时间里，我们一直没有一个很好的数学模型去表示自然语言背后的语义含义。而这节的方法叫 LSI，Latent Semantic Analysis（潜在语义分析）——有时也叫 Latent Semantic Indexing，所以有时叫 LSA、有时叫 LSI，俩是一回事。它描述的就是：怎么把高维空间的文本数据降到低维空间，使语义的信息比较清晰地呈现出来——这是第一次有人用 Semantic（语义）这个词，第一次尝试用词向量的模型来描述语义。

走一遍：从 6 维到 2 维

在这个矩阵里，每一个数据（一行）是文本，每一个维度是一个词——每个特征代表「某个词对于这个文本是否出现、出现的频率」。6 个词就是 6 维。把它做奇异值分解，取最大的两维，从 6 维压到二维。板书上的数字：第一个文档压成了 (−0.888, 0.46)，第二个是 (0.759, 0.652)——每个原先六维的数据，都变成了一个新的二维数据。

（此处原为课件图示，现用文字描述：二维平面上标出 6 个文档点与查询点 Q1。第一、二、三个文档聚在平面一侧——它们讲石头；第四、五、六个聚在另一侧——它们讲音乐。查询 Q1（rock marble）落在石头那一簇里，紧挨着 D1。）

在二维数据上我们就可以在向量空间里把它们画出来：一二三讲的是石头，四五六讲的是音乐。现在 query「rock marble」进来——离 Q1 最近的是 D1，它们在讲关于石头的问题。不同的语义，在这个空间里通过数学的方式比较完美地呈现了出来；至少，我们可以通过它们的余弦夹角来求语义之间的距离。

GT 的 6×6 矩阵完整数字在板书上不可考，下面用一组可手算的缩微数据，把整个流程走一遍（词聚合成两列：石头类 rock/granite/marble 与音乐类 music/song/band；文档 4 篇 + 1 个 query）：

import math

# LSI 缩微例（讲义设定：GT 的 6 词 x 6 文档矩阵完整数字不可考，此处把词聚合为两列频次）
A = [
    [6, 5, 0, 1],    # 石头类词频: D1  D2 | D3  D4
    [0, 1, 6, 5],    # 音乐类词频: D1  D2 | D3  D4
]
docs = ["D1", "D2", "D3", "D4"]
query = [2, 0]       # 新文本 "rock marble": 只有石头类词出现

n_terms, n_docs = len(A), len(A[0])

# 对词空间做 SVD: A*A^T 是 2x2, 求根公式手算特征分解 (上一节同款)
AAt = [[sum(A[i][k] * A[j][k] for k in range(n_docs)) for j in range(n_terms)] for i in range(n_terms)]
print("A * A^T (2x2):")
for row in AAt:
    print(" ", row)

trace = AAt[0][0] + AAt[1][1]
det = AAt[0][0] * AAt[1][1] - AAt[0][1] * AAt[1][0]
disc = math.sqrt(trace * trace - 4 * det)
lambda1, lambda2 = (trace + disc) / 2, (trace - disc) / 2
u1 = [AAt[0][1], lambda1 - AAt[0][0]]
u1 = [x / math.sqrt(u1[0]**2 + u1[1]**2) for x in u1]
u2 = [-u1[1], u1[0]]
print(f"sigma1 = {math.sqrt(lambda1):.4f}, sigma2 = {math.sqrt(lambda2):.4f}   (两个奇异值)")
print(f"concept1 u1 = [{u1[0]:.4f}, {u1[1]:.4f}]   (总词频方向)")
print(f"concept2 u2 = [{u2[0]:.4f}, {u2[1]:.4f}]   (石头-音乐对比方向)")

print("二维概念空间坐标 (concept1, concept2):")
coords = []
for j in range(n_docs):
    d = [A[0][j], A[1][j]]
    c = [d[0]*u1[0] + d[1]*u1[1], d[0]*u2[0] + d[1]*u2[1]]
    coords.append(c)
    print(f"  {docs[j]}: ({c[0]:.4f}, {c[1]:.4f})")
qc = [query[0]*u1[0] + query[1]*u1[1], query[0]*u2[0] + query[1]*u2[1]]
print(f"  query(rock marble): ({qc[0]:.4f}, {qc[1]:.4f})")

def cosine(a, b):
    dot_ab = sum(x*y for x, y in zip(a, b))
    na = math.sqrt(sum(x*x for x in a)); nb = math.sqrt(sum(y*y for y in b))
    return dot_ab / (na * nb)

print("query 与各文档的余弦相似度 (越大越近):")
for j in range(n_docs):
    print(f"  query vs {docs[j]}: {cosine(qc, coords[j]):.4f}")

A * A^T (2x2):
  [62, 10]
  [10, 62]
sigma1 = 8.4853, sigma2 = 7.2111   (两个奇异值)
concept1 u1 = [0.7071, 0.7071]   (总词频方向)
concept2 u2 = [-0.7071, 0.7071]   (石头-音乐对比方向)
二维概念空间坐标 (concept1, concept2):
  D1: (4.2426, -4.2426)
  D2: (4.2426, -2.8284)
  D3: (4.2426, 4.2426)
  D4: (4.2426, 2.8284)
  query(rock marble): (1.4142, -1.4142)
query 与各文档的余弦相似度 (越大越近):
  query vs D1: 1.0000
  query vs D2: 0.9806
  query vs D3: 0.0000
  query vs D4: 0.1961

三个看点：

• query 落进了石头簇：它跟纯石头文档 D1 的余弦恰好是 1.0000（两个向量完全同向——都只含石头类词），跟次近的 D2 有 0.9806；跟纯音乐文档 D3 的余弦是 0.0000——石头与音乐，在这个概念空间里正交，语义上毫无关系。
• 第二个概念轴 u₂ = (−0.7071, 0.7071) 就是「石头减音乐」的对比方向：石头文档坐标为负、音乐文档为正，一目了然。
• 自然语言的问题一般就是高维问题——词典的词数一般是 2 万维、3 万维，多的到 10 万维，所以非常需要映射到低维空间。（图像数据的维度也很高，但自然语言确实是典型的高维问题。）

三种方法的递进

把这个方法放回课程的大图景里。上次课讲 topic model 之前，先讲了 PLSA；今天讲的是 LSA。三者的关系是：

1. LSA：纯是通过一种代数的方法进行降维——就是本节的 SVD；
2. PLSA：用了一部分概率的形式，给它多加了一个随机变量（第 15 章讲过的三硬币结构）；
3. topic model（LDA）：在 PLSA 的基础之上，又假设了一层新的随机变量（Dirichlet 先验）。

从代数到概率再到多层概率——这条线是自然语言语义建模一路走来的脚印。课件里还给了第二个例子：词是 data、inference、retrieval、brain、lung，前四个文档都含 data/retrieval 这类词，是讲 computer science 的；后两个是 brain（大脑）和 lung（肺），讲医学的。SVD 分解之后得到的 $U$、$\Sigma$、$V^T$ 三个矩阵（这里展示的 $\Sigma$ 已经是降完维的），把每个文本乘上 $V^T$ 就映射到二维空间里去。课件对这组矩阵的解读值得记下：$U$ 是文档-概念相似度矩阵、$V$ 是词-概念相似度矩阵、$\sigma_1$ 是「CS 概念」的强度——$V_{1,1}$ 是词 data 在 CS 概念上的权重，$V_{1,2}=0$ 表示 data 与第二个概念（Medical）零相似。这样的代码网上很多，并不难找到。

从线性到非线性：瑞士卷与麻花

讲到这里，请大家注意：以上所有方法都属于线性变换。我们是把高维数据通过 SVD 也好、PCA 也好，映射到低维空间——乘一个矩阵所做的事情，一定是线性变换。那非线性变换该怎么做？仍然是高维空间数据，怎么给它进行变换？这就要接着往下说了。

首先回到上次课（16.2 节）没讲完的问题：新月形状的聚类。那还只是平面里的复杂结构；如果数据在高维空间里，呈现的特性更复杂——比如像一个瑞士卷（Swiss Roll），我们怎么把瑞士卷展开？如果它是一个互相缠绕的更复杂结构，像一根麻花，怎么解开？瑞士卷至少还是一层互相卷，麻花是两根流缠绕在一起的结果。

流形：高维空间里很薄的一小块

这些数据都涉及一个概念：很多数据是存在于高维空间里的流形（Manifold）。什么意思？比如说，不同人脸的不同角度和光照，在几万维的像素空间里，实际只占非常小的一块——这一小块组成一个有内在关系的区域。怎么学到、找到这一小块流形，对整个问题的分析是非常重要的一环。

怎么理解「非常小的一块」？讲师给了一个思想实验。一张 $50 \times 80$ 的图片有 4000 个像素点。现在随机生成一张 $50 \times 80$ 的图——4000 个点里随机取任何的颜色——你会得到什么？基本就是一个雪花噪点，一种不可知、不可识别的东西。你再随机采一张，还是噪点。任何一张图片都是 4000 维空间里的一个点；随便在 4000 维空间里踩一个点，踩出来的图像都是白花花、乱七八糟的。

但别忘了：生活中确实有很多 $50 \times 80$ 的图片是能看出东西的——有的像太阳，有的是场景，有的是杯子、计算机——我们一看就联想到生活里的东西。能识别的这些图，只占 4000 维空间里非常非常小的一部分，而且这一小部分之间是有关系的：一张杯子的图片做些简单的线性变换，相当于这个点在 4000 维空间里从一个地方走到了附近的地方——这些「附近」连起来，就是一张很薄很薄的流形。你可以把流形理解成高维空间里的一个「面」（严格说是子空间，但比如 4000 维空间里的 3000 维子空间，谁也没法想象——一般把它想成三维空间里的一个曲面就够了）。

Isomap：把距离换成图上最短路

那怎么进行非线性的降维？第一个方法叫 Isomap，它解决的正是瑞士卷的问题。瑞士卷最重要的特征是一层一层卷在一起——关键在于 distance measure：距离怎么衡量，它决定了两个点「是否相近」。

拿一维的卷举例（点 1 到 7 依次排在一条卷曲的弧上）：

• 用欧氏距离：点 1 和点 2 的距离是这样一条直线，1 和 3、1 和 4、1 和 7 的距离各是一条直线——卷内两层之间「抄近道」直接穿过去；
• 用 Dijkstra 算法的距离：点 1 到 2 是一段；1 到 3 必须经过 2 再到 3；1 到 5 必须经过 1→2→3→4→5——一个点到另一个点的距离，是在整个图里的最短路径，沿着弧一步步「连加」出来。

非常简单的思路。这样当我们把三维的卷映射到二维的时候——距离不再是欧氏距离，整个卷就相当于被两个小人从两头拉平了。在这种情况下再做聚类，会得到更好的结果。

回来看 16.2 节的新月问题：想衡量新月一端与另一端两点的距离时，实际上是通过一连串最近邻的接力来衡量的——这个点最近的是它，它的最近的是下一个……想求任何两点的距离，就是这条链上距离的连加。于是我们可以根据这种新的距离来进行聚类——什么都不用变：K-means 算法完全不变，还是找中心点、找最近距离，变的只是距离的运算方式（从欧氏距离换成 Dijkstra 的图上最短路径）。用代码把「换距离」这件事亲眼看到（两个月牙、每个点只与最近邻连边、朴素 Dijkstra）：

import math

# 两个月牙（讲义设定：两条圆弧, 模拟 16.2 节 K-means 失灵的新月数据）
def arc_points(cx, cy, r, angles_deg):
    pts = []
    for a in angles_deg:
        rad = math.radians(a)
        pts.append((cx + r * math.cos(rad), cy + r * math.sin(rad)))
    return pts

upper = arc_points(3, 0, 3, [30, 75, 105, 150])     # 上月牙, 索引 0-3
lower = arc_points(7, 0, 3, [210, 255, 285, 330])   # 下月牙, 索引 4-7
points = upper + lower
n = len(points)

def dist(i, j):
    return math.dist(points[i], points[j])

# 近邻图: 每个点与它欧氏距离最近的 1 个点连边 (课件: 欧氏距离阈值或 top-k 近邻)
edges = set()
for i in range(n):
    nearest = sorted(range(n), key=lambda j: dist(i, j))[1]  # [0]是自己, 取最近 1 个
    edges.add((min(i, nearest), max(i, nearest)))

adj = {i: [] for i in range(n)}
for i, j in sorted(edges):
    w = dist(i, j)
    adj[i].append((j, w))
    adj[j].append((i, w))

# 朴素 Dijkstra: 从 s 出发到各点的图上最短路
def dijkstra(s):
    INF = float('inf')
    d = [INF] * n
    d[s] = 0.0
    visited = [False] * n
    for _ in range(n):
        u, best = -1, INF
        for v in range(n):
            if not visited[v] and d[v] < best:
                best, u = d[v], v
        if u == -1:
            break
        visited[u] = True
        for v, w in adj[u]:
            if d[u] + w < d[v]:
                d[v] = d[u] + w
    return d

d3 = dijkstra(3)   # 上月牙左端点
print(f"同月牙两端 (上月牙3号-0号): 欧式 {dist(3, 0):.4f} | 图最短路 {d3[0]:.4f}")
print(f"跨月牙最近对 (0号-4号): 欧式 {dist(0, 4):.4f} (全图最近的跨弧对) | 图最短路 {d3[4]}")
print(f"跨月牙 (3号-4号): 欧式 {dist(3, 4):.4f} | 图最短路 {d3[4]}")

同月牙两端 (上月牙3号-0号): 欧式 5.1962 | 图最短路 6.1451
跨月牙最近对 (0号-4号): 欧式 3.2297 (全图最近的跨弧对) | 图最短路 inf
跨月牙 (3号-4号): 欧式 5.0000 | 图最短路 inf

结果一目了然：

• 欧氏距离最近的跨月牙对（3.2297），在图上根本不可达（inf）——它们之间没有任何「近邻链」能连过去。两个月牙在近邻图上是两个不连通的分量，天然就是两簇。
• 同一条月牙的两端，欧氏直线 5.1962，图上要沿着弧一步步连加出 6.1451——图距离略长，但它是「沿着流形走」的真实距离。

这就是 Isomap 的全部思想：聚类算法一行不改，把距离的度量换成近邻图上的最短路径，新月、瑞士卷就摊开了。（当然，要映射到低维空间里去，仍然可以用这套图距离去做——课件给的形式是 $R^3 \to R^2$ 的 embedding，即嵌入。）

LLE：用邻居定义的点

当然还有一种方式，叫 LLE，Local Linear Embedding（局部线性嵌入）：把高维空间的流形映射到低维空间时，让它们的线性关系得到保留。这个方法稍稍复杂一点，讲师明说：数学原理今天可能讲不完，先把核心思路给大家描述清楚。

第一步（在高维空间学权重）：在高维空间的流形上，对每个点找到它最近的 $k$ 个近邻（比如对红色点找它最近的 4 个邻居）。我们需要训练一个 $W$（权重）来表示它和邻居之间的关系：给定的这个 $x_i$，用周围那几个邻居点 $x_j$ 的线性组合（乘以权重再相加）来重新构造它，学习一个很好的 $W$，使重构的误差最小——也就是优化这个代价函数 $E(W)$。$W$ 满足一个约束：权重加和等于 1。这一步是可以求出 $W$ 的。

第二步（在低维空间保关系）：把 $m$ 维空间嵌入（embedding）到更低维的空间里去时，希望 $W$ 的关系得到保留：原先的 $x_i$ 映射到新空间里是 $y_i$、$x_j$ 映射过去是 $y_j$——在新的空间里，用以前学到的那些 $W$ 值，同样能用周围的几个邻居重新构造出新的 $y$。找到满足这个条件的映射，使整体误差 $E(Y)$ 最小。

跟 PCA 对照一下就看出味道了：PCA 乘一个矩阵的映射，流形上原有的关系就消失了；而 LLE 的映射是在「每个点都由它周围的点定义」这个约束下进行的——约束下的映射，保留了高维空间的流形特性。（讲师说 Local Linear Embedding 每个公式的具体推导，需要下节课再跟进——下一章的 LLE 专节正是兑现这句话的地方。）

本章收线：无监督学习四件事

这一章（也是这两次课）讲的内容稍稍有点杂，讲师最后自己串了一遍：

1. LDA 的延伸：任何一个 topic 可以变成「一个词的分布」的形式，而这个分布可以用另外一个 Dirichlet 分布表示——给 topic 一个更宏观的定义；
2. 明确了 topic、词、文本三者之间的关系；
3. 无监督学习四件事：第一是 k-means（聚类）、第二是 PCA（主成分分析）、第三是 SVD（奇异值分解）、第四是从线性映射到非线性——矩阵分解之外的 Isomap、LLE 这些办法。

其中 LLE 没有给到具体的算法和数学推导——那是接下来章节的任务。

📝 动手练一练

1. 方法递进辨析：LSA、PLSA、topic model（LDA）三者是什么递进关系？各自的「数学工具箱」里最核心的东西分别是什么？（提示：一个纯代数、一个加了概率和一个随机变量、一个又假设了一层随机变量——结合第 15 章的三硬币结构与 Dirichlet 先验。）

2. 余弦手算：概念空间里，文档 D1 的坐标是 $(4.2426, -4.2426)$、query 的坐标是 $(1.4142, -1.4142)$。手算两者的余弦相似度，并解释为什么它恰好是 1。再用同样方法算 query 与 D3 $(4.2426, 4.2426)$ 的余弦，解释那个 0 意味着什么。

本章小结

SVD 把任何矩阵分解成 $U \Sigma V^T$ 三个矩阵的乘积，$\Sigma$ 的对角线上是奇异值——取最大的几个、其余归零，就得到原矩阵的低维逼近 $\hat{A}$。它最经典的应用是 LSI：词-文档矩阵经 SVD 压进概念空间后，一词多义（rock）能靠「跟谁在低维空间里近」来消歧，语义的距离用余弦夹角来量。LSA（纯代数）→ PLSA（概率+随机变量）→ LDA（再加一层先验）是语义建模的三级台阶。而所有「乘矩阵」的办法都是线性变换；面对流形数据，Isomap 把距离换成近邻图上的最短路（K-means 一行不改），LLE 用「邻居线性组合」的约束保住流形结构。

📋 行动清单

• 手算走查：把 LSI 缩微例的 $A A^T$ 特征分解在纸上重推一遍（trace=124、det=3744、两个特征值 72 和 52），核对本节代码的每一个数。
• 复盘三方法：回到第 15 章的 PLSA 与 LDA 讲义，对照本节的「三方法递进」，画出 LSA→PLSA→LDA 的关系图（各自加了什么、各自牺牲了什么）。
• 预习下一章：聚类、局部线性嵌入（LLE）与分布式随机邻域嵌入——本节欠下的 LLE 公式推导，下一章兑现。
• 思考应用：搜索引擎里输入一个有多义性的词（如「苹果」），想一想 LSI 的概念空间会怎么摆放「水果新闻」「手机新闻」和你的 query？

—— 小象教研组