""" ============================================================================= 第 3 讲 · 代码 3：用 numpy 手写两层神经网络 —— 前向、反向、为什么需要非线性 ============================================================================= 对应讲义：lectures/lecture03_neural_networks_backprop.md 第 3、5 节 ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 代码 01 推导了一层网络的梯度，代码 02 用标量计算图自动求导。真实的网络是"成批"处理数据的： 输入是一个 (N, d) 的矩阵，每层都是矩阵乘法。本文件用 numpy 手写一个**向量化**的两层网络： 前向：H = f(X W1 + b1)，logits = H W2 + b2，损失 = 平均交叉熵 反向：dlogits = (softmax - Y) / N → dW2 = Hᵀ dlogits → dH = dlogits W2ᵀ → dZ = dH ∘ f'(Z) → dW1 = Xᵀ dZ 数据是经典的"螺旋"数据集：三类点各自绕成一条螺旋线，彼此缠绕，**不能用直线分开**。 演示 1：线性分类器（softmax 回归）在螺旋数据上只能做到约一半的准确率 演示 2：两层网络的梯度检验 演示 3：两层网络把三条螺旋完美分开；终端里画出决策边界 演示 4：去掉非线性激活函数会怎样？sigmoid、tanh、ReLU 哪个好训练？ ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python 03_numpy_mlp.py 需要 numpy；装了 matplotlib 会把决策边界图保存到 figures/。运行时间约 5~10 秒。 """ from pathlib import Path import numpy as np FIGURE_DIR = Path(__file__).parent / "figures" def print_header(title): """打印一个带分隔线的标题。""" print() print("=" * 70) print(f" {title}") print("=" * 70) # ============================================================================= # 第 1 部分：数据 # ============================================================================= def make_spiral(n_per_class=100, n_classes=3, noise=0.2, seed=0): """ 生成螺旋数据集：每一类的点沿着一条螺旋线分布，各类螺旋互相缠绕。 【思路】 第 k 类的第 i 个点：半径 r 从 0 均匀增加到 1，角度 θ 从 4k 均匀增加到 4(k+1)，再加上噪声。 点的坐标 = (r·sin θ, r·cos θ)。 【参数】 n_per_class : 每类多少个点 n_classes : 类别数 noise : 角度上的高斯噪声标准差（越大越难分） seed : 随机种子 【返回】 (X, y)：X 形状 (N, 2)，y 形状 (N,)，取值 0..n_classes-1 【NLP 联系】 二维数据看得见、画得出，适合建立直觉。NLP 中的输入（词向量拼接）有成百上千维， 但"线性分类器只能画直线（超平面）、神经网络可以画任意形状的边界"这个道理完全一样。 """ rng = np.random.default_rng(seed) X, y = [], [] for k in range(n_classes): r = np.linspace(0.0, 1.0, n_per_class) theta = np.linspace(4 * k, 4 * (k + 1), n_per_class) + rng.normal(scale=noise, size=n_per_class) X.append(np.stack([r * np.sin(theta), r * np.cos(theta)], axis=1)) y.append(np.full(n_per_class, k)) return np.concatenate(X), np.concatenate(y) # ============================================================================= # 第 2 部分：激活函数 # ============================================================================= # 每个激活函数都给出 (f, f')：f' 用前向的结果来表示，避免重复计算。 ACTIVATIONS = { "relu": (lambda z: np.maximum(0, z), lambda z, h: (z > 0).astype(float)), "tanh": (np.tanh, lambda z, h: 1 - h ** 2), "sigmoid": (lambda z: 1 / (1 + np.exp(-z)), lambda z, h: h * (1 - h)), "linear": (lambda z: z, lambda z, h: np.ones_like(z)), # "没有激活函数" } # ============================================================================= # 第 3 部分：两层神经网络 # ============================================================================= class TwoLayerNet: """ 两层全连接网络：输入 → 隐藏层（激活函数 f）→ 输出 logits → softmax。 【参数】（构造函数） n_in, n_hidden, n_out : 输入维度、隐藏单元数、类别数 activation : "relu"、"tanh"、"sigmoid" 或 "linear" seed : 随机种子 【边界】 - 权重用 N(0, 1/n_in) 初始化（方差和输入维度成反比），偏置初始化为 0。 代码 05 会专门讨论初始化为什么重要。 """ def __init__(self, n_in, n_hidden, n_out, activation="relu", seed=0): rng = np.random.default_rng(seed) self.params = { "W1": rng.normal(scale=np.sqrt(1.0 / n_in), size=(n_in, n_hidden)), "b1": np.zeros(n_hidden), "W2": rng.normal(scale=np.sqrt(1.0 / n_hidden), size=(n_hidden, n_out)), "b2": np.zeros(n_out), } self.f, self.df = ACTIVATIONS[activation] def forward(self, X): """ 前向传播，返回 logits，并把反向传播需要的中间结果存进 self.cache。 【参数】X : np.ndarray (N, n_in) 【返回】np.ndarray (N, n_out) """ p = self.params Z = X @ p["W1"] + p["b1"] # (N, n_hidden) H = self.f(Z) logits = H @ p["W2"] + p["b2"] # (N, n_out) self.cache = (X, Z, H) return logits def loss_and_grads(self, X, y, l2=0.0): """ 计算平均交叉熵损失和所有参数的梯度。 【思路】（讲义第 5.5 节，每一步都是代码 01 里验证过的规则） dlogits = (softmax(logits) - Y) / N softmax + 交叉熵 dW2 = Hᵀ dlogits，db2 = dlogits 按行求和 线性层（对批里每个样本的梯度求和） dH = dlogits W2ᵀ 线性层对输入的梯度 dZ = dH ∘ f'(Z) 逐元素激活函数：逐元素相乘 dW1 = Xᵀ dZ，db1 = dZ 按行求和 L2 正则：损失加上 (l2/2)·Σ W²，梯度加上 l2·W（只对权重，不对偏置）。 【参数】 X : np.ndarray (N, n_in) y : np.ndarray[int] (N,) l2 : float，L2 正则系数 【返回】 (loss, grads)：grads 是和 self.params 同结构的字典，每个梯度和对应参数同形状 【Python 小知识】 X.T @ dZ：形状 (n_in, N) @ (N, n_hidden) = (n_in, n_hidden)，和 W1 一样。 这里的矩阵乘法自动完成了"对批里所有样本的梯度求和"。 """ p = self.params logits = self.forward(X) X, Z, H = self.cache N = len(y) shifted = logits - logits.max(axis=1, keepdims=True) probs = np.exp(shifted) probs /= probs.sum(axis=1, keepdims=True) loss = -np.mean(np.log(probs[np.arange(N), y])) loss += 0.5 * l2 * (np.sum(p["W1"] ** 2) + np.sum(p["W2"] ** 2)) dlogits = probs.copy() dlogits[np.arange(N), y] -= 1.0 dlogits /= N grads = {"W2": H.T @ dlogits + l2 * p["W2"], "b2": dlogits.sum(axis=0)} dH = dlogits @ p["W2"].T dZ = dH * self.df(Z, H) grads["W1"] = X.T @ dZ + l2 * p["W1"] grads["b1"] = dZ.sum(axis=0) return float(loss), grads def predict(self, X): """返回每个样本得分最高的类别。""" return np.argmax(self.forward(X), axis=1) def train(model, X, y, lr=1.0, steps=3000, l2=1e-3, log_every=None): """ 全批量梯度下降训练（数据只有几百个点，每步都用全部数据）。 【参数】 model : TwoLayerNet 或 LinearSoftmax X, y : 训练数据 lr : 学习率；steps : 迭代步数；l2 : L2 正则系数 log_every : 每多少步打印一次（None 表示不打印） 【返回】 list[float]，每一步的损失 """ history = [] for step in range(steps): loss, grads = model.loss_and_grads(X, y, l2) for k in model.params: model.params[k] -= lr * grads[k] history.append(loss) if log_every and (step % log_every == 0 or step == steps - 1): acc = np.mean(model.predict(X) == y) print(f" 第 {step:>4d} 步：损失 {loss:.4f}，训练准确率 {acc:.1%}") return history class LinearSoftmax(TwoLayerNet): """ 线性分类器（softmax 回归）：logits = X W + b，没有隐藏层。 【边界】 - 为了复用 train()，它和 TwoLayerNet 有相同的接口（params、loss_and_grads、predict）。 - 它的决策边界只能是直线（在高维中是超平面）。 """ def __init__(self, n_in, n_out, seed=0): rng = np.random.default_rng(seed) self.params = {"W": rng.normal(scale=0.01, size=(n_in, n_out)), "b": np.zeros(n_out)} def forward(self, X): return X @ self.params["W"] + self.params["b"] def loss_and_grads(self, X, y, l2=0.0): logits = self.forward(X) N = len(y) probs = np.exp(logits - logits.max(axis=1, keepdims=True)) probs /= probs.sum(axis=1, keepdims=True) loss = -np.mean(np.log(probs[np.arange(N), y])) + 0.5 * l2 * np.sum(self.params["W"] ** 2) d = probs.copy() d[np.arange(N), y] -= 1.0 d /= N return float(loss), {"W": X.T @ d + l2 * self.params["W"], "b": d.sum(axis=0)} # ============================================================================= # 第 4 部分：画决策边界 # ============================================================================= def ascii_decision_boundary(model, X, y, width=64, height=26): """ 在终端里画出决策边界：背景字符（. - :）表示模型在这个位置会预测哪一类， 数据点用 A、B、C 表示它们的真实类别（预测错的点用小写字母）。 【返回】str """ xs = np.linspace(X[:, 0].min() - 0.1, X[:, 0].max() + 0.1, width) ys = np.linspace(X[:, 1].max() + 0.1, X[:, 1].min() - 0.1, height) grid = np.array([[x, yy] for yy in ys for x in xs]) pred = model.predict(grid).reshape(height, width) canvas = [[".-:"[pred[r, c]] for c in range(width)] for r in range(height)] pred_pts = model.predict(X) for (px, py), label, pl in zip(X, y, pred_pts): c = int(round((px - xs[0]) / (xs[-1] - xs[0]) * (width - 1))) r = int(round((ys[0] - py) / (ys[0] - ys[-1]) * (height - 1))) ch = "ABC"[label] canvas[r][c] = ch if pl == label else ch.lower() return "\n".join(" " + "".join(row) for row in canvas) def save_boundary_png(model, X, y, path, title): """（可选）用 matplotlib 保存决策边界图。没装 matplotlib 时返回 False。""" try: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt except ImportError: return False xx, yy = np.meshgrid(np.linspace(-1.1, 1.1, 300), np.linspace(-1.1, 1.1, 300)) Z = model.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape) fig, ax = plt.subplots(figsize=(5, 5)) ax.contourf(xx, yy, Z, alpha=0.3, levels=[-0.5, 0.5, 1.5, 2.5]) ax.scatter(X[:, 0], X[:, 1], c=y, s=10) ax.set_title(title) fig.tight_layout() fig.savefig(path, dpi=110) plt.close(fig) return True # ============================================================================= # 演示 # ============================================================================= def demo_linear(X, y): """演示 1：线性分类器。""" print_header("演示 1：线性分类器（softmax 回归）在螺旋数据上") model = LinearSoftmax(2, 3) train(model, X, y, lr=1.0, steps=1000, log_every=250) print("\n 决策边界（. - : 是模型预测的区域，A B C 是真实类别，小写 = 预测错误）：") print(ascii_decision_boundary(model, X, y)) print(""" 观察：线性分类器只能用直线把平面切成三块，而三条螺旋互相缠绕，怎么切都会切错一大片。""") return model def demo_gradient_check(): """演示 2：两层网络的梯度检验。""" print_header("演示 2：两层网络的梯度检验") rng = np.random.default_rng(1) Xs, ys = rng.normal(size=(5, 3)), rng.integers(0, 4, size=5) for act in ["relu", "tanh", "sigmoid"]: net = TwoLayerNet(3, 6, 4, activation=act, seed=2) _, grads = net.loss_and_grads(Xs, ys, l2=0.1) worst = 0.0 for k, P in net.params.items(): num = np.zeros_like(P) it = np.nditer(P, flags=["multi_index"]) while not it.finished: i = it.multi_index old = P[i] P[i] = old + 1e-5 lp, _ = net.loss_and_grads(Xs, ys, l2=0.1) P[i] = old - 1e-5 lm, _ = net.loss_and_grads(Xs, ys, l2=0.1) P[i] = old num[i] = (lp - lm) / 2e-5 it.iternext() rel = np.max(np.abs(num - grads[k])) / max(np.max(np.abs(num)), np.max(np.abs(grads[k])), 1e-12) worst = max(worst, rel) print(f" 激活函数 {act:<8s}：所有参数梯度的最大相对误差 {worst:.1e}") assert worst < 1e-5 print(" （ReLU 在 0 点不可导，数值梯度恰好跨过 0 时会有误差；随机数据下几乎不会发生。）") def demo_mlp(X, y): """演示 3：两层网络。""" print_header("演示 3：两层神经网络（100 个 ReLU 隐藏单元）") model = TwoLayerNet(2, 100, 3, activation="relu") train(model, X, y, lr=1.0, steps=3000, log_every=500) print("\n 决策边界：") print(ascii_decision_boundary(model, X, y)) FIGURE_DIR.mkdir(exist_ok=True) if save_boundary_png(model, X, y, FIGURE_DIR / "lecture03_spiral_mlp.png", "2-layer ReLU network"): print(" （已保存图片：figures/lecture03_spiral_mlp.png）") print(""" 观察：隐藏层把每个点映射到一个 100 维的新空间，在新空间里三类点变得可以用"直线"（超平面）分开； 映射回原来的二维平面，边界就成了弯曲的螺旋形。神经网络 = 自动学习一个让问题变简单的特征表示。""") def demo_activations(X, y): """演示 4：不同的激活函数。""" print_header("演示 4：激活函数的影响（隐藏层 100 个单元，学习率 1.0，同样训练 3000 步）") print() for act in ["linear", "sigmoid", "tanh", "relu"]: model = TwoLayerNet(2, 100, 3, activation=act) hist = train(model, X, y, lr=1.0, steps=3000) acc = np.mean(model.predict(X) == y) print(f" {act:<8s}：第 300 步损失 {hist[299]:.3f}，最终损失 {hist[-1]:.3f}，训练准确率 {acc:.1%}") print(""" 观察： - linear（没有激活函数）：两个线性层叠起来还是线性函数 W2(W1 x) = (W2 W1) x，和演示 1 的线性分类器一样差。 **非线性激活函数是神经网络能力的来源。** - sigmoid 学得最慢：它的导数最大只有 0.25，而且输出不以 0 为中心； tanh 的导数最大是 1；ReLU 在正半轴导数恒为 1，梯度传得最顺畅，也最便宜（只要比较一下大小）。 这是 ReLU 成为默认选择的主要原因（讲义第 3.4 节）。""") if __name__ == "__main__": X, y = make_spiral() print(f"螺旋数据集：{len(X)} 个点，{len(set(y))} 类，每个点 2 维") demo_linear(X, y) demo_gradient_check() demo_mlp(X, y) demo_activations(X, y)