Ch 2 · 线性代数扫盲(NumPy 实战)

一句话:1-2 小时过完向量、矩阵、分解这些符号,让你见到时认得出、大致懂,要求能默写每一行证明。

§1 目标

Ch 1 之后我们已经会写 Python。现在翻到 Ch 3 那篇 PyTorch 教程,代码段里全是 model = nn.Linear(in, out)x @ W.T + bloss.backward(),公式那一栏全是 W x + bsoftmax(qk^T / √d) · vtorch.linalg.svd(...)。这些符号背后其实就是本章要见一遍的矩阵运算——我们这章的目标就是把这一层符号认得,让它们不再是天书

具体说,跑完本章,我们达到三件事。其一——见到 PyTorch / Tensorflow / sklearn 文档里跳出来的 A @ BA.Tnp.linalg.solvetorch.linalg.svd@ 等符号时,认得出它们在做哪类运算(矩阵乘法、转置、解方程组、SVD),能顺着函数名把「公式」跟「代码」对上号。其二——见到支持向量机(SVM)、主成分分析(PCA)、梯度下降这些后续章节要上的概念时,知道它们背后跑的就是本章要见的工具:SVM 涉及二次规划(线性方程组的扩展),PCA 涉及特征分解或 SVD,梯度下降每一步都涉及梯度——梯度是导数,机器学习里几乎所有导数都是矩阵-向量运算其三——为 Ch 3 的最优化(梯度下降 / 牛顿法 / KKT)、Ch 4 的 SVM、Ch 5 的 PCA 打一层最基本的铺垫,后面遇到具体公式不会「符号卡住」。

承诺的深度对齐到「见过 / 认得 / 扫过」——承诺让你像数学系那样流利地手写证明,承诺让你口算一个 5 × 5 矩阵的特征值。CNN 课的本质是「会用库解决问题」,不是「数值线性代数」。

不承诺:

  • ❌ 矩阵证明每一步都自己推——那是另一门课的活儿。
  • ❌ 见到任何奇异矩阵都能瞬间判定为什么不可逆——靠经验,本章只给直觉。
  • ❌ 口算 5 × 5 矩阵的特征值——下面会反复强调「调库就行」。

§2 为什么先做这一章

Ch 1 之后,你会用 dataclass 写训练配置、用 pathlib 读 JSON、用 uv 装包。这些都是「会写代码」的一面。

但 CNN 的本质是矩阵运算 + 梯度下降。一旦翻到 PyTorch 教程,公式那栏全是 WbW^T x、那个 softmax(...),demo 代码那栏全是「一行定义 layer、一行算 loss、一行 backward」。你不认得符号,就连 demo 代码是「在算什么」都不知道

Ch 2 的策略因此克制:只够后续章节用,不抢「线性代数课」的活儿。三个 demo、十几个函数、几十个测试用例,跑完得了「认得」深度。这些不是一遍就能吃透的事——Ch 2 过一遍,Ch 3 + 在用到的位置加深。

本章不教什么(明确范围):

  • 矩阵理论的严格证明(行列式怎么算、特征值为什么存在)——CNN 课不做证明,只做直觉。
  • 稀疏矩阵、QR 分解、张量运算——后续章节不会集中用到,本章不引入。
  • 复数域、向量空间公理化——超出 CNN 课范围。
  • 自己手写一份数值稳定的 SVD / 矩阵分解实现——见到库调 np.linalg.svd 即可。

§3 概念

接下来 8 个子节,3.1 介绍本章要跑的三个 demo,3.2 把 ndarray 基础过一遍,3.3-3.7 是五类运算(向量 / 矩阵 / 解方程组 / 特征分解 / SVD),3.8 把范数统一下。

3.1 本章的项目

本章要跑的是三个纯 numpy 的短 demo,加起来不到 200 行:

  1. demo 1 — 向量与范数 (code.ex1_vectors_and_norms):向量长度(L1 / L2)、内积、线性组合、正交判定。
  2. demo 2 — 矩阵运算 (code.ex2_matrix_ops):@ 乘法、转置、行列式、求逆(奇异矩阵的行为)。
  3. demo 3 — 矩阵分解 (code.ex3_decompositions):解线性方程组、特征分解、SVD 重建(低秩近似)。

每个 demo 暴露四五个函数 + 一个 main() 入口。每个函数都有对应单元测试,共 57 个测试用例(详见 5 节)。三块跑完,我们过了「把数学公式翻译成 numpy」这套最小子集的一关。

3.2 numpy ndarray 基础

Ch 1 没碰过 numpy,先把三件事过一遍——下面所有运算都建立在 ndarray 之上。

ndarray 是 numpy 的核心数据结构,与 Python list 的根本区别:所有元素同一个数据类型、按多维数组组织。深度学习里,一个 batch 的图像就是一个四维 ndarray(数量、高度、宽度、通道),整张梯度表就是个二维 ndarray。两个 ndarray 相加 / 相乘,按逐元素规则广播。

三件必须记住的事:

.shape 是维度的元组。形状 (3, 4) 的矩阵有 3 行 4 列、12 个元素。形状是后续所有运算的前提——大多数「莫名 ValueError」都来自形状不匹配。

.dtype 是元素类型。np.array([1, 2, 3]) 默认 int64,np.array([1.0, 2.0]) 默认 float64。SVD、解线性方程组都要 float——传 int 数组会立即报错或返回垃圾。

Broadcasting——形状兼容的小数组自动扩展到与大数组相同的形状:

import numpy as np

x = np.array([1.0, 2.0, 3.0])
print(x + 1.0)              # [2. 3. 4.] (标量被广播成 [1. 1. 1])
print(np.ones((3, 3)) + x)  # 每行都加 x,得到 shape (3, 3)

广播规则:从右往左对齐维度,某一维不等就要么是 1 要么不存在,否则报错。3.4 节矩阵减向量就是这条规则的直接应用。

3.3 向量:点积、范数、线性组合、正交

向量是「一组有序数」的几何对象。在 CNN 里,一张图像「展平」(flatten) 后就是一个长度为像素数的长向量;神经网络一层的输出也是一个向量。CNN 里绝大多数运算都建立在两个向量之间的关系上——这一节先见这层关系。

点积 dot(v, u) 是最基础的关系——两个向量的「投影乘积」,公式上等价于逐元素乘再求和。直观上:方向越接近,点积越大;正交时点积为 0;反向时点积为负

v = np.array([3.0, 4.0])
u = np.array([4.0, -3.0])
print(np.dot(v, u))    # 0.0  (正交)
print(np.dot(v, v))    # 25.0 (跟自己的点积 = L2 范数的平方)

L2 范数是「向量的长度」,定义为「v 跟 v 的点积开根号」。一般化的 p-范数有 L1(绝对值之和)、L∞(最大元素绝对值)等,统称 p-范数。numpy 用一个函数覆盖:

v = np.array([3.0, 4.0])
print(np.linalg.norm(v, ord=2))  # 5.0 (3-4-5 三角形)
print(np.linalg.norm(v, ord=1))  # 7.0 (|3| + |4|)

线性组合 c1 * v1 + c2 * v2 + ... 是「用一组基向量构造新向量」的通用操作。n 维空间里任意向量都能写成 n 个基向量的线性组合——这就是「坐标」的代数含义。

e1 = np.array([1.0, 0.0])
e2 = np.array([0.0, 1.0])
print(2.0 * e1 + 3.0 * e2)   # [2. 3.] ← 标准基下的坐标

正交是两个向量夹角为 90°,等价于「点积的绝对值在数值精度内 ≈ 0」。CNN 里正交初始化的权重矩阵训得更快——这是后话,本章先认得这词。

3.4 矩阵:乘法、转置、求逆、行列式

矩阵是「一组有序向量」或「一次线性变换」的几何对象。CNN 里每一层权重 W 就是一个矩阵;整张图像就是一个矩阵(高 × 宽,或 通道 × 高 × 宽)。理解矩阵 = 理解神经网络。

矩阵乘法要求左矩阵的列数 = 右矩阵的行数,结果的第 (i, j) 元素是「左矩阵第 i 行点积右矩阵第 j 列」。「输入向量 x 左乘 W 等于线性变换后的输出」在 numpy 里就是 @:

W = np.array([[4.0, 7.0], [2.0, 6.0]])
I = np.eye(2)
print(W @ I)   # [[4. 7.] [2. 6.]] (单位矩阵是矩阵乘法里的「1」)

转置 .T 把行变成列、列变成行,形状 (M, N) 变成 (N, M),沿对角线「翻折」。线性回归公式 W^T x 里那个 T 就是转置。

W = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])
print(W.T.shape)   # (2, 3)

行列式 det(W) 是方阵的一个标量。直观含义:矩阵代表的「线性变换」对体积的拉伸倍数。det(W) == 0 时,矩阵是奇异的——降维、不可逆、列向量平行。奇异矩阵没有逆

逆矩阵 inv(W) 是满足 W @ inv(W) = I 的矩阵。只有方阵 + 非奇异时才有逆。生产里几乎从不显式求逆——求 W 的逆 乘 b 直接用 np.linalg.solve(W, b)(见 3.5 节),数值上稳定得多。

W = np.array([[4.0, 7.0], [2.0, 6.0]])
print(np.linalg.det(W))       # 10.0 (非零 → 可逆)
print(np.linalg.inv(W))       # [[0.6, -0.7], [-0.2, 0.4]]
print(W @ np.linalg.inv(W))   # ≈ [[1, 0], [0, 1]] (单位阵 + 浮点尾巴)

3.5 解线性方程组 W x = b

很多问题能列成「已知系数矩阵 W、目标向量 b、求变量向量 x」的形式。例如:三道工序分别用时 x1, x2, x3 分钟,产量总和是 b 道——3 未知数 3 方程。numpy 里就一行 np.linalg.solve(W, b)

几何上,每一行方程定义 n 维空间里一个超平面,多个平面交于一点就是解。W 可逆(行列式非零)时,解存在且唯一;否则要么无解、要么无穷多解(对应奇异矩阵)。

W = np.array([[3.0, 2.0], [1.0, 2.0]])
b = np.array([7.0, 5.0])
x = np.linalg.solve(W, b)
print(x)        # [1. 2.]
print(W @ x)    # [7. 5.] ← 代回原方程

注意 np.linalg.solve(W, b) 不是 np.linalg.inv(W) @ b,虽然数学上等价。前者数值稳定得多——直接求逆再相乘,在矩阵接近奇异时会放大误差。

3.6 特征分解:特征值与特征向量

对某些方阵 W,存在特殊方向 x(特征向量)使得 W @ xx 只在长度上有差别——矩阵把这些方向「拉长 / 缩短」但不旋转。这个方向叫特征向量,拉伸倍数叫特征值:

W @ x == 特征值 * x

特征值大的方向是矩阵「最活跃」的方向。PCA(主成分分析)就是找数据协方差矩阵的最大几个特征值对应的特征向量——这些方向是数据方差最大的几个「轴」。

W = np.array([[2.0, 1.0], [1.0, 2.0]])   # 对称矩阵
eigvals, eigvecs = np.linalg.eig(W)
print(eigvals)   # [3.+0.j 1.+0.j]  ← 实对称 → 实特征值,以 complex 形式返回
print(eigvecs)   # 特征向量按列排列

坑点 1:np.linalg.eig任何方阵都返回 complex128 类型数组,即使特征值是实数——对称矩阵永远有实特征值,但 numpy 还是给 complex。比较前要 .real 取实部:

eigvals, _ = np.linalg.eig(W)
print(eigvals.dtype)        # complex128
print(eigvals.real)         # [3. 1.]
print(np.sort(eigvals.real))  # 对实部排序(特征值的顺序依赖实现)

坑点 2:只有方阵有特征分解。对称 → 实特征值 + 正交特征向量;非对称 → 可能出现复特征值。本章只过到「认得」层级,严格证明留在线性代数课里。

3.7 SVD:任意矩阵都能分解

特征分解只能用在方阵;任意 m × n 矩阵都能做奇异值分解 (SVD):

A == U @ 对角矩阵(对角线上排奇异值) @ V 转置

即拆成三块:第一个 m × k 的「左奇异向量」(按列)、第二个 k × k 对角矩阵(对角线上排奇异值,降序)、第三个 k × n 的「右奇异向量转置」。其中 k = min(m, n),每个奇异值都 ≥ 0。

SVD 的杀手锏是低秩近似:取前 j 个奇异值(j < k)重建出来的矩阵,是这个矩阵的所有 rank-j 近似里 Frobenius 范数误差最小的(Eckart-Young 定理)。这是 PCA、推荐系统、图像压缩背后的核心数学工具:

import numpy as np
M = np.random.default_rng(seed=0).standard_normal((4, 4))
U, S, Vh = np.linalg.svd(M, full_matrices=False)
print(S)            # 4 个奇异值,降序
approx = U[:, :2] @ np.diag(S[:2]) @ Vh[:2, :]
print(np.linalg.norm(M - approx, ord="fro"))   # Frobenius 重建误差

注意 np.linalg.svd 的默认值在历史上不同版本不一致(True / False 都见过),full_matrices=True 会让 U 是 m×m、Vh 是 n×N,被零行零列「撑胖」,既浪费内存也容易让切片索引混乱。默认就显式传 full_matrices=False,UVh 都只有 k 列。这是写 SVD 代码的默认姿势。

3.8 范数:L1 / L2 / Frobenius

范数是「把一个向量 / 矩阵量成一个长度」的非负函数。最常见的三类:

  • L1 范数(向量):绝对值之和。np.linalg.norm(v, ord=1)
  • L2 范数(向量):欧几里得长度。np.linalg.norm(v, ord=2)
  • Frobenius 范数(矩阵):所有元素的平方和开根号,等于「把矩阵展平成向量后的 L2 范数」。np.linalg.norm(A, ord="fro")

CNN 里常见的「L2 正则化」就是给权重矩阵加一个「W 的 L2 范数平方」形式的惩罚项(加进 loss 里);Ch 3 的最优化章节也会大量出现这些范数。

v = np.array([3.0, 4.0])
print(np.linalg.norm(v, ord=1))         # 7.0
print(np.linalg.norm(v, ord=2))         # 5.0
A = np.array([[1.0, 2.0], [3.0, 4.0]])
print(np.linalg.norm(A, ord="fro"))     # sqrt(1+4+9+16) ≈ 5.477

§4 动手

三个 demo 的 main() 跑出来。先确认环境:uv sync 在本章目录跑通即可——Ch 1 已经把 uv 讲清楚了。

4.1 把环境搭起来

cd articles/02-linear-algebra-numpy
uv sync

uv sync 不报错,说明依赖齐了。下面的命令都在这个目录下执行。

4.2 demo 1 — 向量与范数

跑:

$ python -m code.ex1_vectors_and_norms
L2 norm of [3. 4.] = 5.0000
dot(v, u) = 0.0 (理论上 0)
orthogonal(v, u) = True
linear_combine([2.0, 3.0], bases) = [2. 3.]

逐行解读:

  • L2 norm of [3. 4.] = 5.0000 — 3-4-5 三角形的斜边长,验证范数定义「3 的平方 + 4 的平方开根号」。
  • dot(v, u) = 0.0v = [3, 4]u = [4, -3] 正交(3*4 + 4*(-3) = 0),所以内积为 0。
  • orthogonal(v, u) = True — 上面那一点用「绝对值 < 容差(默认 1e-8)」判定正交。
  • linear_combine([2.0, 3.0], bases) = [2. 3.][2, 3] 当成标准基下的坐标还原回向量,得到 [2, 3](因为标准基就是 (1, 0)(0, 1))。

四个输出里 dot(v, u) = 0.0精确等于 0,不是「数值上接近 0」——这是挑了 [3, 4][4, -3] 这对手工对出来的向量,直接相乘抵消。如果换成 [1, 0] vs [0, 1],结果同样是 0.0,但实际是 1*0 + 0*1 = 0,肉眼可见的正交。

4.3 demo 2 — 矩阵运算

跑:

$ python -m code.ex2_matrix_ops
matmul(A, I) = [[4.0, 7.0], [2.0, 6.0]]
det(A) = 10.0000
inv(A) = [[0.6, -0.7000000000000001], [-0.2, 0.4]]
verify A @ inv(A) = [[0.9999999999999998, -1.1102230246251565e-16], [-1.1102230246251565e-16, 1.0]]

逐行解读:

  • matmul(A, I) = [[4.0, 7.0], [2.0, 6.0]] — 单位矩阵是矩阵乘法的「1」,A @ I == A
  • det(A) = 10.0000A = [[4, 7], [2, 6]],行列式 4*6 - 7*2 == 10,非零 → A 可逆。
  • inv(A) = [[0.6, -0.7000000000000001], [-0.2, 0.4]] — 伴随矩阵 / 行列式。注意 -0.7000000000000001 不是错,是浮点表示 0.7 的精度尾巴——下面验算会回到单位阵。
  • verify A @ inv(A) = [[0.99..., -1.11e-16], [-1.11e-16, 1.0]] — 对角线 ≈ 1,非对角线是「极小量」(1e-16 数量级),数值上就是单位阵。「是不是单位阵」不该用 == 判等,应该用 assert_allclose(atol=1e-10)

4.4 demo 3 — 矩阵分解

跑:

$ python -m code.ex3_decompositions
解 A @ x = b: x = [1. 2.]
验证 A @ x = [7. 5.]
SVD rank-2 重建误差 ||A - approx||_F = 1.1986

逐行解读:

  • 解 A @ x = b: x = [1. 2.] — 解 [[3, 2], [1, 2]] @ x = [7, 5],得到 x = [1, 2]
  • 验证 A @ x = [7. 5.] — 代回原方程,左边确实等于右边。
  • SVD rank-2 重建误差 ||A - approx||_F = 1.1986 — 一个 4×4 随机矩阵用前 2 个奇异值重建后,Frobenius 范数误差是 1.1986。这是信息损失的下界——丢掉一半奇异值后的最小可能误差(Eckart-Young)。

§5 验证

三个 demo 都能跑 + 单元测试全绿。完整命令清单:

# 1. 环境健康
uv sync

# 2. 三个 demo 入口(预期输出贴在 §4 各小节)
python -m code.ex1_vectors_and_norms
python -m code.ex2_matrix_ops
python -m code.ex3_decompositions

# 3. 单元测试
python -m pytest tests/ -v

预期输出:

  • 三个 demo 的输出已在 §4 贴出。
  • pytest 报告:57 passed(57 = 三个 demo 全部测试方法数,含 @pytest.mark.parametrize 展开后的总数;非展开大约 30 个测试方法)。

如果 pytest 失败,看测试方法名直接定位段:TestMatmul::* 失败 → matmul 函数出问题;TestSvdReconstruct::* 失败 → svd_reconstruct 函数出问题,以此类推。

§6 回顾

本章带你扫过(仅「见过」深度):

  • numpy ndarray 基础(.shape.dtype、broadcasting)
  • 向量:点积、L1 / L2 范数、线性组合、正交判定(含容差)
  • 矩阵:@.Tnp.linalg.detnp.linalg.inv、奇异矩阵行为
  • 线性方程组 W x = b:np.linalg.solve 的形态与「不要先求逆」的告诫
  • 特征分解(特征值 × 特征向量形式):np.linalg.eig、对称 vs 非对称、complex128 习惯
  • SVD:三块矩阵相乘(np.linalg.svd(full_matrices=False))、低秩重建、Eckart-Young 直觉
  • 范数回顾:L1 / L2 / Frobenius
  • 三个 demo 共 12 个函数,57 个测试用例(uv syncpython -m pytest tests/ -v 全绿)

本章不展开(后续章节用到时再加深):

  • 特征值的数值稳定性(幂迭代、QR 迭代等算法)
  • 稀疏矩阵(scipy.sparse)——CNN 课用 dense 输入,暂不引入
  • QR 分解、Cholesky、LUP——某些数值算法里替代 SVD 的工具,本章用不到
  • SVD 的严格推导和奇异值的几何含义
  • 复数域的特征分解(非对称矩阵的复特征值、特征多项式)
  • 范数的对偶性
  • 伪逆 np.linalg.pinv——非方阵最小二乘解的标准工具
  • 广播的全部规则(右对齐、维度 1 扩展、不匹配报错)
  • batched np.einsum——CNN 中常用的多维张量乘法简写,Ch 3 以后才上

§7 下篇预告

下一章会用本章的 numpy 工具 + PyTorch 的 autograd,首次进入最优化。我们会见到:

  • 梯度下降的最小例子:L 对 θ 求偏导,沿负梯度方向走一步。
  • 为什么 CNN 用 SGD / Adam 而不是牛顿法——海森矩阵在百万参数下没法算。
  • 一个能在笔记本上跑通的最小训练循环——届时见到 loss.backward(),能认出「这是 autograd 在算梯度」。

不会推论文,不会讲理论推导——只讲「跑起来需要知道什么」。