1
2
3
4
5
6
7
8

规格化数的表示式（尾数隐含前导 1，故实际精度是"尾数位 + 1"）：
  x = (−1)^s × 1.f₂ × 2^(e − bias)         f₂ 是尾数位构成的二进制小数
  有效位数 p = 尾数位 + 1 = 24（FP32）/ 53（FP64）
  相邻可表示数之间的最大相对间距 = 2^(−p)，称为单位舍入误差 u：
      FP32 u = 2^(−24) ≈ 5.96e−8       FP64 u = 2^(−53) ≈ 1.11e−16
"机器 epsilon" 的两种常见口径（读论文前先确认作者用哪一种）：
  口径一（1 与下一个更大浮点数之差）：ε = 2^(−(p−1))，FP32 ≈ 1.19e−7、FP64 ≈ 2.22e−16
  Python 的 sys.float_info.epsilon 用的是口径一，值为 2.22e−16

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

0.1 × 2 = 0.2 → 取整 0      0.1₁₀ = 0.0001100110011001100…₂
0.2 × 2 = 0.4 → 取整 0              └─ 0011 无限循环 ─┘
0.8 × 2 = 1.6 → 取整 1       规格化后：1.1001100110011…₂ × 2^(−4)
1.6 × 2 = 1.2 → 取整 1
0.2 × 2 = 0.4 → 取整 0   ⟲ 回到第二行，循环开始

FP64 中实际存储的值：0.1 → 0.1000000000000000055511151231257827…（偏大）
                     0.2 → 0.2000000000000000111022302462515654…（偏大）
                     0.3 → 0.2999999999999999888977697537484346…（偏小）
  0.1 + 0.2 = 0.3000000000000000444089209850062616… ≠ 0.3

两个误差没有抵消：0.1 与 0.2 的存储值都偏大，相加后误差叠加，
恰好落在 0.3 的存储值的另一侧。只有分母是 2 的幂的十进制小数
（0.5、0.25、0.125）才能在二进制下精确表示。

1
2
3

设 x̂ = x(1 + δ₁)、ŷ = y(1 + δ₂)，|δ_i| ≤ u，则 x̂ − ŷ = (x − y) + (xδ₁ − yδ₂)，
绝对误差 ≤ u·(|x| + |y|)，相对误差 ≤ u·(|x| + |y|) / |x − y|
⟹ 当 x ≈ y 时分母 → 0，相对误差被放大到无法接受，这就是灾难性相消。

1
2
3

        e^(x_i)              e^(−m)·e^(x_i)            e^(x_i − m)
  ─────────────────  =  ─────────────────────────  =  ─────────────────
    Σ_j e^(x_j)          e^(−m)·Σ_j e^(x_j)          Σ_j e^(x_j − m)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

记 μ = E[X]，σ² = Var(X)，则 E[X²] = μ² + σ²，单遍法要算 (μ² + σ²) − μ²：
两个操作数量级都是 μ²，真正的结果却只有 σ²。
  放大因子 = (μ² + σ²)/σ² = 1 + (μ/σ)²      ⟹ 相对误差 ~ u·(1 + (μ/σ)²)

本例取 x = {1e8+1, 1e8+2, 1e8+3}，真值 Var = 2/3 ≈ 0.667。
单遍法得 E[X²] ≈ 1e16 + 4e8 + 4.667 与 (E[X])² ≈ 1e16 + 4e8 + 4，
两者量级都是 1e16，而 FP64 在该量级的相邻间距 ulp = 2^(53−52) = 2
（1e16 ≈ 2^53.15，二进制指数为 53），故两个操作数各自被舍入到最近的偶数，
差值只能取 0 或 2，根本得不到 0.667。本例 μ/σ = 1.2e8，放大因子 ~1.5e16，
乘 u ≈ 1.1e−16 得相对误差 ~1.6 > 1：结果完全失真，"方差"甚至可能为负。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

def kahan_sum(values):
    """补偿求和：把每次加法丢掉的低位收集到 comp 里，下一轮再补回去"""
    total = 0.0
    comp = 0.0                       # 上一轮被舍入丢掉的部分
    for x in values:
        y = x - comp                 # 先把欠账补上
        t = total + y                # 这一步可能丢掉低位
        # (t - total) 在浮点下是精确的（两个相近数相减无误差，见 5.2.2），
        # 因此它能无损地还原 total + y 时被丢掉的量；减去 y 就得到误差本身
        comp = (t - total) - y
        total = t
    return total

1
2
3
4
5
6
7
8
9

唯一性：设 p(x)、q(x) 都是次数 ≤ n 且过全部 n+1 个节点的多项式，
  令 r(x) = p(x) − q(x)。则 deg r ≤ n，且对每个节点
      r(x_i) = p(x_i) − q(x_i) = y_i − y_i = 0
  即 r 有 n+1 个互不相同的根。而次数不超过 n 的非零多项式最多有 n 个根，
  所以只能是 r ≡ 0，即 p = q。

存在性（构造）：定义 Lagrange 基函数 ℓ_i(x) = Π_{j ≠ i}(x − x_j)/(x_i − x_j)。
  由构造，ℓ_i(x_i) = 1（分子分母相同），ℓ_i(x_j) = 0（j ≠ i 时分子含 (x_j − x_j)）。
  于是 L(x) = Σ_{i=0}^{n} y_i·ℓ_i(x) 在每个节点取值为 y_i，存在性得证。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

N(x) = f[x₀] + f[x₀,x₁](x − x₀) + f[x₀,x₁,x₂](x − x₀)(x − x₁) + …
                … + f[x₀,…,x_n](x − x₀)(x − x₁)…(x − x_{n−1})

差商递推：f[x_i] = y_i，f[x_i,…,x_{i+k}] = (f[x_{i+1},…,x_{i+k}] − f[x_i,…,x_{i+k−1}])/(x_{i+k} − x_i)
差商表（自上而下、自左而右地算，左下对角线就是 Newton 系数）：
      x_i    f[x_i]        一阶                   二阶
      0      1
                           (2−1)/(1−0) = 1
      1      2                                   (3−1)/(2−0) = 1
                           (5−2)/(2−1) = 3
      2      5

验证（真题 【2026-08-05 选择 Q9】的三个点 (0,1), (1,2), (2,5)）：
  N(x) = 1 + 1·(x − 0) + 1·(x − 0)(x − 1) = 1 + x + x² − x = x² + 1
  N(1.5) = 2.25 + 1 = 3.25  ✓    与直接解方程组 a = 1, b = 0, c = 1 的答案一致

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

Jacobi：左端只留 D，整轮全用上一轮的旧值
  D x^(k+1) = b − (L + U)x^(k) ⟹ x^(k+1) = −D⁻¹(L + U)x^(k) + D⁻¹b
  迭代矩阵 B_J = −D⁻¹(L + U)；分量形式 x_i^(k+1) = (b_i − Σ_{j≠i} a_ij x_j^(k))/a_ii

Gauss-Seidel：把 L 也挪到左端，算新分量时立刻用上本轮已更新的值
  (D + L)x^(k+1) = b − U x^(k) ⟹ x^(k+1) = −(D + L)⁻¹U x^(k) + (D + L)⁻¹b
  迭代矩阵 B_GS = −(D + L)⁻¹U
  分量形式 x_i^(k+1) = (b_i − Σ_{j<i} a_ij x_j^(k+1) − Σ_{j>i} a_ij x_j^(k))/a_ii
                       └─ 用本轮新值 ─┘          └─ 用上轮旧值 ─┘

SOR（逐次超松弛，松弛因子 ω）：在 Gauss-Seidel 的更新上做外推
  (D + ωL)x^(k+1) = [(1 − ω)D − ωU]x^(k) + ωb ⟹ B_ω = (D + ωL)⁻¹[(1 − ω)D − ωU]
  ω = 1 时退化为 Gauss-Seidel；0 < ω < 2 是收敛的必要条件。

1
2
3
4

B_J = −D⁻¹(L + U) 的对角元为 0，非对角元为 −a_ij/a_ii，故
  ‖B_J‖_∞ = max_i Σ_{j≠i} |a_ij| / |a_ii| < 1        （正是严格对角占优的定义）
而谱半径不超过任何诱导范数：ρ(B_J) ≤ ‖B_J‖_∞ < 1，故迭代收敛。
Gauss-Seidel 的证明更复杂，但结论相同：严格对角占优时它也收敛。

1
2
3

设 Ax = b，而实际算了 A(x + δx) = b + δb。两式相减得 A·δx = δb，即 δx = A⁻¹δb，
故 ‖δx‖ ≤ ‖A⁻¹‖·‖δb‖；又由 b = Ax 得 ‖b‖ ≤ ‖A‖·‖x‖，即 1/‖x‖ ≤ ‖A‖/‖b‖。
两式相乘：‖δx‖/‖x‖ ≤ ‖A‖·‖A⁻¹‖·‖δb‖/‖b‖ = cond(A)·‖δb‖/‖b‖

1
2
3
4
5
6
7
8

min_x ‖Ax − b‖² + λ‖x‖²
对 x 求导置零：2Aᵀ(Ax − b) + 2λx = 0 ⟹ (AᵀA + λI)x = Aᵀb ⟹ x = (AᵀA + λI)⁻¹Aᵀb
为什么能改善条件数：AᵀA 的奇异值是 σ_i²，而 AᵀA + λI 的是 σ_i² + λ，于是
  cond(AᵀA + λI) = (σ_max² + λ)/(σ_min² + λ) < σ_max²/σ_min² = cond(AᵀA)
λ 越大条件数越接近 1，但偏差也越大——这就是第 4 章 4.6.3 节的偏差—方差权衡。
第 8 章会从贝叶斯角度再解释：λ‖x‖² 等价于对权重施加均值为 0、方差为 1/(2λ) 的
高斯先验，岭回归就是该先验下的 MAP 估计。同一件事在数值分析里叫"改善条件数"，
在统计学习里叫"正则化"，在贝叶斯里叫"先验"，三者是同一个数学对象。

1
2
3
4
5

① 绝对容差 |a − b| ≤ EPS，适用于量级已知且相差不大的场景（坐标、评分）
② 相对容差 |a − b| ≤ EPS·max(1, |a|, |b|)，适用于量级未知或跨度很大的场景
   （max 中的 1 是为了兼容 a、b 接近 0）
③ 组合判据（Python math.isclose 的默认行为）：|a − b| ≤ max(REL·max(|a|,|b|), ABS)
选择 EPS 的唯一原则：EPS 必须【大于】累加噪声（约 n·u·尺度），同时远【小于】真实需要区分的最小间隙。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

EPS = 1e-9   # 比累加噪声（约 1e-15）大一个量级，又远小于真实可分辨的代价间隙

def pick_best(costs):
    """按 (代价, 序号) 字典序取最优：代价相等时保留序号靠前的那个"""
    best_idx, best_cost = -1, 0.0
    for i, cost in enumerate(costs):
        # 判据是 cost < best_cost − EPS 而不是 cost < best_cost：
        # 只有"超出容差地更优"才改判，数学上相等的两者由先遍历到的（序号小的）胜出
        if best_idx < 0 or cost < best_cost - EPS:
            best_idx, best_cost = i, cost
    return best_idx