""" ============================================================================= CS349E · 前置讲 P3 · 03_rounding_and_accumulation.py —— 舍入方式与累加精度：低精度矩阵乘法的误差从哪来 ============================================================================= 【这个脚本要回答的问题】 一次舍入的误差很小（半个 ulp），但一个矩阵乘法要做几千次加法，误差会不会越积越多？ 本脚本用 numfmt 的"每次加法都舍入"的模拟器回答： 演示 1：四种舍入方式（就近取偶、截断、四舍五入、随机舍入）——单次误差和**系统性偏差**。 演示 2：停滞（stagnation）：在 BF16 里一直加 1，加到 256 就加不动了；随机舍入怎样救回来。 演示 3：FP8 矩阵乘法的累加器：输入是 FP8，乘积累加到 FP32 / "只有 14 位"/ BF16 / FP16 里，K 越大误差越大； DeepSeek-V3 的"每 128 个元素提升到 FP32"能把误差降多少。 演示 4：求和算法：顺序、两两（pairwise）、Kahan 补偿——同样在 BF16 / FP16 里，误差差几个数量级。 【怎么读输出】全部是确定性的模拟（固定种子），每次运行结果相同。 "14 位累加器"是一个**玩具模型**：DeepSeek-V3 报告 H800 上 FP8 GEMM 的累加"只保留约 14 位"， 但没有说明具体的舍入方式，所以这里同时给出"就近舍入"和"截断"两个版本，真实硬件的行为不在本脚本的判断范围内。 【运行】python3 03_rounding_and_accumulation.py （约 3 秒） python3 03_rounding_and_accumulation.py --quick （约 1 秒，演示 3 只算 K ≤ 4096） 【依赖】numpy；从 ../lecture01 导入 console；同目录的 numfmt。 """ import os import sys import numpy as np HERE = os.path.dirname(os.path.abspath(__file__)) sys.path.insert(0, os.path.join(HERE, "..", "lecture01")) from console import print_header, print_observations, print_table # noqa: E402 from numfmt import BF16, FP4_E2M1, FP8_E4M3, FP16, FP32, FloatFormat, accumulate # noqa: E402 QUICK = "--quick" in sys.argv # 玩具模型："只保留约 14 位有效数字"的累加器 = FP32 的指数范围 + 13 位尾数（加上隐含的 1 共 14 位） ACC14 = FloatFormat("ACC14", 8, 13, 127, "ieee") # --------------------------------------------------------------------------- # 演示 1：舍入方式 # --------------------------------------------------------------------------- def demo1_rounding_modes(n=200_000, seed=0): """ 第一部分：把几个数舍入到 FP4 E2M1（间隔大，看得清楚），比较四种方式。随机舍入给出 1 万次的平均值。 第二部分：把 [1, 2) 里均匀分布的 n 个数舍入到 BF16（ulp = 2^−7），统计误差的**平均值**（偏差）和均方根， 单位都是 ulp。 【为什么偏差比误差大小更重要】一次舍入误差 ≤ 半个 ulp，单独看很小。但如果每次都往同一个方向偏（截断总是变小）， 几千次累加后误差就线性增长；无偏的误差会互相抵消，只按 √n 增长。 【随机舍入 stochastic rounding】以"x 离下面那个数的距离 / ulp"为概率向上取。单次误差更大（最多一个 ulp）， 但期望值正好是 x——无偏。低精度训练（Gupta et al. 2015）和一些 AI 芯片的权重更新用它。 """ print_header("演示 1：四种舍入方式——单次误差与系统性偏差") modes = [("nearest_even", "就近，平局取偶（IEEE 默认）"), ("toward_zero", "截断（向 0）"), ("nearest_away", "就近，平局远离 0（四舍五入）"), ("stochastic", "随机舍入")] xs = [0.25, 0.75, 1.25, 2.5, 3.5, 5.0, -2.5, 2.4] rng = np.random.default_rng(seed) rows = [] for key, label in modes: cells = [label] for x in xs: if key == "stochastic": r = FP4_E2M1.round(np.full(10_000, x), rounding="stochastic", rng=rng) cells.append(f"均值 {r.mean():.3f}") else: cells.append(f"{float(FP4_E2M1.round(x, rounding=key)):g}") rows.append(cells) print_table(["舍入到 FP4 E2M1"] + [f"{x:g}" for x in xs], rows) x = rng.uniform(1.0, 2.0, n) ulp = BF16.eps # [1, 2) 区间里 BF16 的 ulp = 2^-7 rows = [] for key, label in modes: err = (BF16.round(x, rounding=key, rng=rng) - x) / ulp rows.append([label, f"{err.mean():+.4f}", f"{np.sqrt(np.mean(err ** 2)):.3f}", f"{np.abs(err).max():.3f}"]) print(f"\n[1, 2) 里均匀分布的 {n:,} 个数舍入到 BF16（单位：ulp = 2^−7）：") print_table(["方式", "平均误差（偏差）", "均方根误差", "最大误差"], rows) print_observations([ "平局（正好在两个数中间）时各方式才有区别：0.25 在 0 和 0.5 中间，就近取偶给 0，四舍五入给 0.5；" "2.5 → 2（取偶）或 3（远离 0）。截断永远往 0 的方向走（2.4 → 2，5 → 4）。", "截断的平均误差是 −0.5 ulp：每次都偏小半个 ulp，这是系统性偏差，累加时会线性积累（演示 3）。", "就近取偶的偏差约为 0：平局时一半向上、一半向下。均匀分布的数平局极少，所以四舍五入和取偶在这里几乎一样，" "但在整数、低精度格式里平局很常见（例如 FP4），这时'远离 0'会产生偏差。", "随机舍入的均方根误差（约 0.41 ulp）比就近舍入（约 0.29 ulp）大，但它是无偏的，而且对'比半个 ulp 还小的'更新也有反应" "（演示 2）。", ]) # --------------------------------------------------------------------------- # 演示 2：停滞 # --------------------------------------------------------------------------- def demo2_stagnation(trials=200, seed=0): """ 在各种格式里从 0 开始一直加 1（或一直加 0.01），看累加器停在哪里。 【为什么会停】累加器 s 越来越大，它的 ulp 也越来越大。当 1 < ulp(s) / 2 时，s + 1 舍入回 s：加了等于没加。 BF16 有 8 位有效数字：到 256 时 ulp = 2，256 + 1 = 257 正好是平局，取偶回到 256。 【随机舍入】s + 1 以 1/ulp 的概率进位成 s + ulp，期望仍是 s + 1。这里重复 trials 次，报告平均值和标准差。 """ print_header("演示 2：停滞——在低精度累加器里一直加一个小数") rows, stop = [], {} for fmt, inc, n in ((BF16, 1.0, 1000), (BF16, 0.01, 10_000), (FP16, 1.0, 10_000), (FP16, 0.01, 10_000)): vals = np.full((1, n), float(fmt.round(inc))) s_rne = accumulate(vals, fmt)[0] stop[(fmt.name, inc)] = s_rne sr = accumulate(np.repeat(vals, trials, axis=0), fmt, rounding="stochastic", seed=seed) exact = vals.sum() rows.append([fmt.name, f"{n:,} 次 + {inc:g}", f"{exact:g}", f"{s_rne:g}", f"{sr.mean():.1f} ± {sr.std():.1f}"]) print_table(["累加器", "操作", "精确答案", "就近舍入", f"随机舍入（{trials} 次的均值 ± 标准差）"], rows) print(f"\n对照：FP32 累加器要加到 2^24 = {2 ** 24:,} 才会停滞（ulp 变成 2）。") print_observations([ f"BF16 里加 1 加到 {stop[('BF16', 1.0)]:g} 就停了（8 位有效数字：256 的 ulp 是 2，256 + 1 是平局、取偶回到 256），" f"加 0.01 停在 {stop[('BF16', 0.01)]:g}。FP16 有 11 位有效数字，加 1 停在 {stop[('FP16', 1.0)]:g}，" f"加 0.01 停在 {stop[('FP16', 0.01)]:g}——都是's 大到增量不足半个 ulp'的地方。", "这类'把很多小数加到一个大数上'的运算无处不在：点积的累加、求均值、softmax 的分母、优化器更新权重（w += lr × g）。" "所以累加器的精度必须比输入高——混合精度的核心原则：低精度存、低精度乘，高精度加。", "随机舍入把停滞救回来了：平均值接近精确答案（代价是有随机波动，标准差见表）。BF16 权重更新用随机舍入，" "就是为了让小于半个 ulp 的更新'按概率'生效。", ]) # --------------------------------------------------------------------------- # 演示 3：FP8 矩阵乘法的累加器 # --------------------------------------------------------------------------- def demo3_fp8_gemm_accumulation(T=256, seed=0): """ 模拟 FP8 矩阵乘法的一个输出元素 c = Σ_k a_k b_k：a、b 是 FP8 E4M3（标准正态随机数量化而来）， 乘积在 FP32 里是精确的（两个 4 位有效数字的乘积只有 8 位），误差全部来自累加。 【累加器】 FP32：标准做法； ACC14 就近 / 截断：玩具模型（见文件开头）； ACC14 截断 + 每 128 提升：每累加 128 个乘积，把部分和加到一个 FP32 累加器里（DeepSeek-V3 的 N_C = 128）； BF16 / FP16：16 位累加器（有的推理库为了省寄存器会这样做）。 【指标】对 T 个输出：范数相对误差 ||ĉ − c|| / ||c||，以及逐元素相对误差 |ĉ − c| / |c| 的中位数。 不用"最大相对误差"：c 恰好接近 0 的元素，相对误差可以任意大，最大值主要由这种元素决定。 """ Ks = (512, 4096) if QUICK else (512, 4096, 16384) print_header(f"演示 3：FP8 矩阵乘法的累加器——输入 FP8 E4M3，误差全部来自累加（每个 K 模拟 {T} 个输出）") configs = [("FP32", FP32, {}), ("ACC14，就近舍入", ACC14, {}), ("ACC14，截断", ACC14, {"rounding": "toward_zero"}), ("ACC14，截断 + 每 128 提升到 FP32", ACC14, {"rounding": "toward_zero", "chunk": 128}), ("ACC14，就近 + 每 128 提升到 FP32", ACC14, {"chunk": 128}), ("BF16", BF16, {}), ("FP16", FP16, {})] rng = np.random.default_rng(seed) results = {} for K in Ks: a = FP8_E4M3.round(rng.standard_normal((T, K))) b = FP8_E4M3.round(rng.standard_normal((T, K))) p = a * b exact = p.sum(axis=1) for name, acc, kw in configs: s = accumulate(p, acc, **kw) norm_rel = np.linalg.norm(s - exact) / np.linalg.norm(exact) med_rel = np.median(np.abs(s - exact) / np.abs(exact)) results[(name, K)] = (norm_rel, med_rel) rows = [] for name, _, _ in configs: cells = [name] for K in Ks: nr, mr = results[(name, K)] cells.append(f"{nr * 100:.3g}% / {mr * 100:.2g}%") rows.append(cells) print_table(["累加器"] + [f"K = {K}：范数 / 中位数" for K in Ks], rows) r_t = results[("ACC14，截断", 4096)][0] r_tp = results[("ACC14，截断 + 每 128 提升到 FP32", 4096)][0] r_n = results[("ACC14，就近舍入", 4096)][0] r_np = results[("ACC14，就近 + 每 128 提升到 FP32", 4096)][0] print_observations([ "FP32 累加几乎精确：误差主要只来自输入的 FP8 量化（这里已经把它当作'精确答案'的一部分）。", f"ACC14（就近）在 K = 4096 时误差 {r_n * 100:.3g}%，截断版 {r_t * 100:.3g}%：截断的偏差每次都朝同一方向，" "随 K 线性积累；就近舍入的误差互相抵消，增长慢得多。DeepSeek-V3 报告 H800 上 K = 4096 时" "'最大相对误差接近 2%'（没有说明误差的定义和舍入方式），和两个玩具模型给出的量级范围一致：" "'只保留约 14 位'足以让长累加的误差变得不可忽略。", f"每 128 个元素提升到 FP32 后，K = 4096 的误差降到 {r_tp * 100:.3g}%（截断）/ {r_np * 100:.3g}%（就近），" "而且几乎不再随 K 增长：低精度累加器里最多只有 128 个数，误差被'重置'。这就是 DeepSeek-V3 的做法，" "代价是每 4 次 WGMMA 要把部分和搬到 CUDA core 上加一次（他们用两个 warpgroup 交替来掩盖这个开销）。", "BF16 累加器在 K = 4096 时误差到了百分之几，K 更大时继续变坏：16 位累加器只适合很短的累加。", ]) # --------------------------------------------------------------------------- # 演示 4：求和算法 # --------------------------------------------------------------------------- def demo4_summation_algorithms(n=10_000, seed=0): """ 同样的 n 个数（[0, 1) 里均匀分布的正数——和一直在变大，最容易暴露顺序求和的问题）， 在 BF16 / FP16 累加器里用三种算法求和。 【三种算法的误差界】（u = 单位舍入误差 = eps / 2；Higham, *Accuracy and Stability of Numerical Algorithms*） 顺序求和：相对误差 ≲ (n − 1)·u —— 随 n 线性增长； 两两求和：相对误差 ≲ ⌈log2 n⌉·u —— GPU 的树形归约、分块求和天然接近这个； Kahan 求和：相对误差 ≲ 2u + O(n·u²) —— 几乎和 n 无关，代价是每个数多 3 次加减法。 这些是**最坏情况**的界；随机数据的误差通常小得多（误差有正有负，按 √n 增长）。 """ print_header(f"演示 4：求和算法——{n:,} 个 [0, 1) 里的正数，在 16 位累加器里求和") rng = np.random.default_rng(seed) rows = [] for fmt in (BF16, FP16): x = fmt.round(rng.uniform(0.0, 1.0, n))[None, :] exact = x.sum() best = float(fmt.round(exact)) # 累加器格式能表示的、离精确答案最近的数 u = fmt.eps / 2 for method, bound in (("sequential", (n - 1) * u), ("pairwise", np.ceil(np.log2(n)) * u), ("kahan", 2 * u)): s = accumulate(x, fmt, method=method)[0] rows.append([fmt.name, method, f"{s:.6g}", f"{exact:.6g}", f"{best:g}", f"{abs(s - exact) / exact:.2e}", f"{bound:.1e}"]) s32 = accumulate(x, FP32)[0] rows.append(["（对照）FP32", "sequential", f"{s32:.6g}", f"{exact:.6g}", "—", f"{abs(s32 - exact) / exact:.2e}", f"{(n - 1) * FP32.eps / 2:.1e}"]) print_table(["累加器", "算法", "结果", "精确答案", "累加器格式里最近的数", "相对误差", "最坏情况的界"], rows) print_observations([ "顺序求和在 BF16 里完全失败（停在 256：之后加的数都不到半个 ulp），在 FP16 里停在 2048，误差都是几十个百分点。", "两两求和把误差降了三个数量级：每个数只参与 log2(n) ≈ 14 次加法，而且相加的两个数大小相近。" "BF16 的两两求和结果 4992 已经是 BF16 能表示的、离精确答案最近的数（5000 附近 BF16 的间隔是 32）。" "GPU 上的并行归约（P2 §6）天然就是两两求和——并行不仅更快，而且更准。", "Kahan 求和在 FP16 里给出 5060，正好是'精确答案舍入到 FP16'的结果（比两两求和的 5056 更接近精确答案）：" "补偿变量 c 记住了每次被舍掉的低位。它的代价是 4 倍的加法，在 memory-bound 的内核里几乎免费" "（多出来的运算藏在访存时间里）。", "实测误差远小于最坏情况的界（顺序求和的界甚至大于 1，没有意义）——界是给'每一次舍入都朝同一方向'的最坏情况准备的。", ]) if __name__ == "__main__": demo1_rounding_modes() demo2_stagnation() demo3_fp8_gemm_accumulation() demo4_summation_algorithms()