附录A：数值稳定性速查表

数值稳定性是大规模矩阵计算的基石。即使算法在理论上正确，数值误差的累积也可能导致计算结果完全失真。本附录汇总了保证数值稳定性的关键技术和实用准则，为实际应用提供快速参考。特别地，我们关注那些在现代机器学习和科学计算中频繁出现但容易被忽视的稳定性问题。

1. 条件数与误差传播

1.1 条件数的多种视角

矩阵条件数 $\kappa(\mathbf{A})$ 衡量了问题的敏感度：

  • 谱条件数：$\kappa_2(\mathbf{A}) = |\mathbf{A}|_2 |\mathbf{A}^{-1}|_2 = \frac{\sigma_{\max}}{\sigma_{\min}}$
  • Frobenius条件数：$\kappa_F(\mathbf{A}) = |\mathbf{A}|_F |\mathbf{A}^{-1}|_F$
  • 分量条件数：考虑输入扰动的结构性

实用估计方法：

  • 使用 condest 类函数进行快速估计（基于1-范数）
  • 通过幂迭代估计最大奇异值
  • 利用Lanczos迭代同时估计最大和最小奇异值

1.2 误差传播的黄金法则

对于线性系统 $\mathbf{A}\mathbf{x} = \mathbf{b}$：

$$\frac{|\Delta \mathbf{x}|}{|\mathbf{x}|} \leq \kappa(\mathbf{A}) \left( \frac{|\Delta \mathbf{A}|}{|\mathbf{A}|} + \frac{|\Delta \mathbf{b}|}{|\mathbf{b}|} \right)$$ 经验准则：

  • 若 $\kappa(\mathbf{A}) \approx 10^k$，则期望损失 $k$ 位有效数字
  • 在双精度下，$\kappa > 10^{16}$ 意味着结果可能完全不可靠
  • 对于迭代方法，条件数还影响收敛速度

1.3 结构化误差分析

前向误差 vs 后向误差：

  • 前向误差：$|\hat{\mathbf{x}} - \mathbf{x}|$ （计算结果与真实解的差距）
  • 后向误差：找到 $\Delta \mathbf{A}$，使得 $(\mathbf{A} + \Delta \mathbf{A})\hat{\mathbf{x}} = \mathbf{b}$

后向稳定算法的优势：

  • 计算结果是某个微扰问题的精确解
  • 若原问题良态，则后向稳定保证前向误差小
  • 大多数LAPACK算法都是后向稳定的

2. 稳定算法速查

2.1 线性系统求解

直接法稳定性排序（从最稳定到最不稳定）：

  1. 带选主元的LU分解 (GEPP) - 增长因子通常为 $O(n^{0.5})$ - 实践中非常稳定

  2. QR分解 - 无条件后向稳定 - 计算成本约为LU的2倍

  3. Cholesky分解（仅适用于正定矩阵） - 不需要选主元 - 自动检测非正定性

  4. 不带选主元的LU分解 - 可能数值不稳定 - 仅在特殊结构下使用

迭代法的稳定性考虑：

  • CG法：条件数影响收敛速度，但每步都保持数值稳定
  • GMRES：理论上有限步收敛，但舍入误差可能导致停滞
  • BiCGSTAB：可能出现崩溃，需要重启策略

2.2 特征值计算

标准算法的稳定性：

  1. 对称矩阵： - Jacobi方法：最稳定，能计算高精度小特征值 - QR算法（带Wilkinson位移）：标准选择 - 分治算法：并行性好，稳定性略低

  2. 非对称矩阵： - Schur分解 + 特征值提取：最稳定 - 直接QR迭代：可能不收敛 - 幂方法：仅适用于主特征值

病态特征值问题的处理：

  • 使用多精度算术
  • 考虑伪谱分析
  • 对于接近的特征值，计算不变子空间而非单个特征向量

2.3 矩阵分解的稳定性

SVD计算：

  • Golub-Reinsch算法：标准方法，双对角化+QR迭代
  • Jacobi SVD：更精确但更慢
  • 随机化SVD：大规模问题的首选，但需要过采样

低秩近似的数值考虑：

  • 截断SVD：最优但计算昂贵
  • 随机化方法：使用过采样参数 $p \geq 5$
  • CUR分解：保持稀疏性，但数值性质较差

2.4 矩阵函数的稳定计算

矩阵指数 $e^{\mathbf{A}}$：

  • Padé近似 + scaling and squaring：通用方法
  • 对于正规矩阵，使用特征分解
  • Krylov方法适用于大规模稀疏情况

矩阵平方根 $\mathbf{A}^{1/2}$：

  • Newton迭代：需要好的初值
  • Schur方法：稳定但计算密集
  • 对于正定矩阵，使用Cholesky分解

矩阵对数 $\log(\mathbf{A})$：

  • 反scaling and squaring
  • 需要 $|\mathbf{A} - \mathbf{I}| < 1$
  • Padé近似的阶数选择很关键

3. 常见数值陷阱与补救措施

3.1 病态问题的识别

预警信号：

  • 残差小但解的扰动大
  • 迭代方法收敛极慢或停滞
  • 不同算法给出显著不同的结果

诊断工具：

  • 计算条件数估计
  • 检查矩阵的奇异值分布
  • 使用扰动分析验证结果

3.2 溢出/下溢的预防

缩放策略：

  • 行平衡：使每行的范数相近
  • 列平衡：使每列的范数相近
  • 对称缩放：保持对称性的同时改善条件数

具体技巧：

行缩放：D₁AD₂，选择D₁和D₂使得缩放后矩阵元素量级均衡
范数计算：使用稳定的实现（如BLAS中的_nrm2）
避免显式形成A^T A

3.3 舍入误差的控制

求和的稳定算法：

  • Kahan求和：补偿舍入误差
  • 成对求和：递归地两两相加
  • 排序求和：从小到大累加

内积计算：

  • 使用BLAS的_dot函数（可能使用扩展精度）
  • 对于高精度需求，使用多精度算术库
  • 考虑误差自由变换（Error-Free Transformations）

3.4 迭代算法的停止准则

相对残差准则： $$\frac{|\mathbf{b} - \mathbf{A}\mathbf{x}_k|}{|\mathbf{b}|} < \text{tol}$$ 加权残差准则（考虑条件数）： $$\frac{|\mathbf{b} - \mathbf{A}\mathbf{x}_k|}{|\mathbf{A}| |\mathbf{x}_k| + |\mathbf{b}|} < \text{tol}$$ 实用建议：

  • 设置最大迭代次数防止无限循环
  • 监控残差的下降历史
  • 对于病态问题，使用更保守的停止准则

4. 特殊场景的稳定性考虑

4.1 稀疏矩阵计算

填充(fill-in)与稳定性的权衡：

  • 最小度排序：减少填充但可能影响稳定性
  • 近似最小度（AMD）：平衡填充和数值性质
  • 嵌套剖分：保证理论界但实现复杂

不完全分解的稳定性：

  • ILU(0)：可能不稳定，需要对角修正
  • ILUT：带阈值，更稳定但填充更多
  • Modified ILU：保持某些数学性质

4.2 分布式计算中的误差累积

通信诱导的误差：

  • 异步更新带来的不一致性
  • 部分和的累积误差
  • 不同处理器的舍入差异

缓解策略：

  • 使用确定性的归约操作
  • 定期同步和误差修正
  • Allreduce操作的稳定实现

4.3 混合精度计算

迭代精化(Iterative Refinement)：

  1. 低精度计算近似解
  2. 高精度计算残差
  3. 低精度求解修正方程
  4. 更新解并重复

精度选择准则：

  • 工作精度：单精度(FP32)或半精度(FP16)
  • 残差计算：至少双精度(FP64)
  • 累加器：考虑使用扩展精度

4.4 GPU计算的数值特性

与CPU的差异：

  • 不同的舍入模式
  • 融合乘加(FMA)的使用
  • 并行归约的非确定性

最佳实践：

  • 使用cuBLAS/cuSolver的稳定例程
  • 注意warp级别的同步
  • 对于高精度需求，考虑使用张量核心的混合精度模式

5. 本章小结

数值稳定性是大规模矩阵计算成功的关键。本附录总结了以下核心概念：

  1. 条件数是理解数值敏感性的基础工具，它量化了输入扰动对输出的放大效应
  2. 后向稳定性是评估算法质量的金标准，确保计算结果是某个邻近问题的精确解
  3. 算法选择需要权衡稳定性、效率和问题结构，没有万能的最优算法
  4. 误差控制需要综合考虑舍入、截断和累积效应，特别是在迭代和分布式场景中
  5. 特殊硬件（如GPU）和混合精度计算带来新的稳定性挑战，需要专门的处理策略

关键公式汇总：

  • 条件数定义：$\kappa(\mathbf{A}) = |\mathbf{A}| |\mathbf{A}^{-1}|$
  • 误差传播：$\text{相对误差} \leq \kappa \times \text{相对扰动}$
  • 机器精度：双精度 $\epsilon \approx 2.2 \times 10^{-16}$
  • 迭代精化收敛条件：$\kappa(\mathbf{A}) \cdot \epsilon < 1$

6. 练习题

基础题

练习 19.1：条件数估计 给定矩阵 $\mathbf{A} = \begin{bmatrix} 1 & 1 \ 1 & 1+\epsilon \end{bmatrix}$，其中 $\epsilon > 0$ 很小。 a) 计算 $\kappa_2(\mathbf{A})$ 作为 $\epsilon$ 的函数 b) 当 $\epsilon = 10^{-8}$ 时，求解 $\mathbf{A}\mathbf{x} = [2, 2+\epsilon]^T$ 预期损失多少位有效数字？

提示：先计算特征值，注意当 $\epsilon$ 很小时的近似。

答案

a) 矩阵 $\mathbf{A}$ 的特征值为：

  • $\lambda_1 = 2 + \epsilon/2 + O(\epsilon^2)$
  • $\lambda_2 = \epsilon/2 + O(\epsilon^2)$

因此：$\kappa_2(\mathbf{A}) = \frac{\lambda_1}{\lambda_2} \approx \frac{4}{\epsilon}$

b) 当 $\epsilon = 10^{-8}$ 时，$\kappa_2(\mathbf{A}) \approx 4 \times 10^8$

预期损失约 8-9 位有效数字。实际解为 $\mathbf{x} = [1, 1]^T$，但数值解可能有显著误差。

练习 19.2：后向误差分析 使用高斯消元法（无选主元）求解： $$\begin{bmatrix} \epsilon & 1 \ 1 & 1 \end{bmatrix} \begin{bmatrix} x_1 \ x_2 \end{bmatrix} = \begin{bmatrix} 1 \ 2 \end{bmatrix}$$ 其中 $\epsilon = 10^{-20}$。分析在双精度下的数值稳定性问题。

提示：追踪消元过程中的数值，注意何时出现灾难性抵消。

答案

消元过程：

  1. 第一步：$m_{21} = 1/\epsilon = 10^{20}$ （溢出风险）
  2. 更新：$a_{22}^{(1)} = 1 - 10^{20}$ （灾难性抵消）
  3. 在浮点运算中：$a_{22}^{(1)} \approx -10^{20}$

回代得到的数值解严重偏离真实解 $x_1 = 1, x_2 = 1$。

使用部分选主元可以避免这个问题：交换两行后，消元因子变为 $\epsilon$，计算稳定。

练习 19.3：迭代精化 考虑线性系统 $\mathbf{A}\mathbf{x} = \mathbf{b}$，其中 $\kappa(\mathbf{A}) = 10^6$。 a) 如果使用单精度（$\epsilon_s \approx 10^{-7}$）求解，迭代精化能否改善结果？ b) 需要多少次迭代才能达到双精度的精度水平？

提示：迭代精化的收敛因子约为 $\kappa(\mathbf{A}) \cdot \epsilon_{\text{work}}$。

答案

a) 收敛因子 $\rho = \kappa(\mathbf{A}) \cdot \epsilon_s = 10^6 \times 10^{-7} = 0.1 < 1$

因此迭代精化会收敛，能够改善结果。

b) 每次迭代误差减少因子约为 0.1。要从单精度误差（$10^{-7}$）改善到双精度（$10^{-16}$），需要误差减少 $10^9$ 倍。

所需迭代次数：$k \geq \frac{\log(10^9)}{\log(10)} = 9$ 次

实践中通常 2-3 次迭代就能获得显著改善。

挑战题

练习 19.4：结构化条件数 考虑三对角矩阵： $$\mathbf{T}_n = \begin{bmatrix} 2 & -1 & & & \ -1 & 2 & -1 & & \ & \ddots & \ddots & \ddots & \ & & -1 & 2 & -1 \ & & & -1 & 2 \end{bmatrix}_{n \times n}$$ a) 证明 $\kappa_2(\mathbf{T}_n) = O(n^2)$ b) 设计一个预条件子使得 $\kappa(\mathbf{M}^{-1}\mathbf{T}_n) = O(1)$ c) 这个结果对有限差分离散化有什么启示？

提示：利用 $\mathbf{T}_n$ 的特征值显式公式，考虑离散二阶导数算子。

答案

a) $\mathbf{T}_n$ 的特征值为：$\lambda_k = 2 - 2\cos\left(\frac{k\pi}{n+1}\right), k = 1, ..., n$

最小特征值：$\lambda_1 \approx \frac{\pi^2}{(n+1)^2}$ 最大特征值：$\lambda_n \approx 4$

因此：$\kappa_2(\mathbf{T}_n) = \frac{\lambda_n}{\lambda_1} \approx \frac{4(n+1)^2}{\pi^2} = O(n^2)$

b) 可以使用如下预条件子：

  • 多重网格方法的粗网格算子
  • 不完全Cholesky分解
  • 循环约化得到的近似逆

这些方法都能使条件数与 $n$ 无关。

c) 启示：

  • 网格加密导致条件数快速增长
  • 简单迭代法（如Jacobi）收敛变慢
  • 多尺度方法是必要的
  • 这解释了为什么PDE求解需要专门的预条件技术

练习 19.5：混合精度策略设计 为以下场景设计混合精度计算策略：

  • 求解大规模稠密线性系统 $\mathbf{A}\mathbf{x} = \mathbf{b}$，$n = 10000$
  • $\mathbf{A}$ 存储需要 800MB（双精度）
  • GPU显存限制为 4GB
  • 需要相对误差 $< 10^{-10}$

要求：最大化性能同时保证数值精度。

提示：考虑分块算法、迭代精化和选择性使用高精度。

答案

优化的混合精度策略：

  1. 存储方案： - $\mathbf{A}$ 存储为FP32（400MB） - 工作空间预留用于分块 - 关键中间量（如主元）保持FP64

  2. 计算流程： ```

  3. FP32: LU分解（带部分选主元）

  4. FP32: 前向/回代得到初始解 x₀
  5. FP64: 计算残差 r = b - Ax₀
  6. FP32: 求解修正方程 Az = r
  7. FP64: 更新 x₁ = x₀ + z

  8. 重复步骤3-5直到收敛 ```

  9. 分块优化： - 块大小选择：平衡显存和计算效率 - 对角块用FP64，非对角块用FP32 - 使用面板分解减少通信

  10. 收敛保证： - 监控相对残差（FP64计算） - 如果停滞，切换到FP64重新分解 - 典型情况下2-3次迭代即可达到目标精度

  11. 性能估计： - 相比纯FP64加速约1.8-2.2倍 - 内存使用减少约45% - 数值精度满足要求

练习 19.6：分布式计算的数值稳定性 在分布式环境中计算 $|\mathbf{A}|_F^2 = \sum_{i,j} a_{ij}^2$，其中矩阵 $\mathbf{A} \in \mathbb{R}^{n \times n}$ 按行分块存储在 $p$ 个处理器上。

a) 分析朴素并行求和的误差累积 b) 设计一个数值稳定的分布式算法 c) 如何处理各处理器上的数据规模严重不平衡的情况？

提示：考虑树形归约和误差补偿技术。

答案

a) 朴素方法的误差分析：

  • 每个处理器局部误差：$O(\epsilon n^2/p)$
  • 归约过程误差：$O(\epsilon \log p)$
  • 最坏情况总误差：$O(\epsilon n^2)$
  • 当部分和数量级差异大时，可能出现灾难性抵消

b) 稳定的分布式算法： ```

  1. 局部计算（每个处理器）：

    • 使用Kahan求和或成对求和
    • 计算局部和 S_i 及误差补偿项 C_i

  2. 全局归约：

    • 使用确定性树形归约
    • 传输(S_i, C_i)对而非单一值
    • 在归约时继续使用补偿求和

  3. 可选的两遍算法：

    • 第一遍：计算规模估计
    • 第二遍：使用缩放避免溢出 ```

c) 负载不平衡的处理：

  • 动态缩放：根据局部最大值选择缩放因子
  • 分层归约：先在相似规模的处理器间归约
  • 使用高精度累加器：如使用128位浮点数
  • 自适应算法： if (max_local_sum / min_local_sum > 10^6) 使用排序归约 else 使用标准树形归约

关键研究方向：

  • 通信与精度的最优权衡
  • 异步算法的误差界
  • 硬件感知的归约策略

练习 19.7：创新算法的稳定性分析 最近提出的"随机递归最小二乘"(Randomized Recursive Least Squares, RRLS)算法通过随机采样和递归更新来求解超定系统。算法核心是： $$\mathbf{H}_{k+1} = \mathbf{H}_k - \frac{\mathbf{H}_k \mathbf{a}_k \mathbf{a}_k^T \mathbf{H}_k}{\mathbf{a}_k^T \mathbf{H}_k \mathbf{a}_k + \lambda}$$ 其中 $\mathbf{H}_k$ 逼近 $(\mathbf{A}^T\mathbf{A})^{-1}$，$\mathbf{a}_k$ 是随机采样的行。

分析该算法的数值稳定性，特别是： a) 何时会出现数值崩溃？ b) 如何改进算法以提高稳定性？ c) 与经典方法相比的优劣？

提示：这是Sherman-Morrison公式的应用，考虑分母接近零的情况。

答案

a) 数值崩溃的条件：

  1. 当 $\mathbf{a}_k^T \mathbf{H}_k \mathbf{a}_k + \lambda \approx 0$
  2. $\mathbf{H}_k$ 失去正定性
  3. 舍入误差累积导致非对称性
  4. 条件数增长失控

b) 稳定性改进方案：

  1. 正定性维护： 每隔m步：H = (H + H^T)/2 如果最小特征值 < ε：H = H + εI

  2. 自适应正则化： λ_k = max(λ_min, τ·trace(H_k)/n) 其中τ ~ 10^(-12)

  3. 数值安全更新： 计算 d = a_k^T H_k a_k + λ if (d < ε·‖H_k‖·‖a_k‖²) 跳过此次更新 else 执行Sherman-Morrison更新

  4. 周期性重启： if (cond(H_k) > threshold) 使用最近的样本重新初始化

c) 与经典方法比较：

优势：

  • 在线/增量更新能力
  • 内存需求 O(n²) vs O(mn)
  • 适合流数据处理

劣势：

  • 收敛速度依赖于采样
  • 数值稳定性较QR分解差
  • 需要精心的参数调节

适用场景：

  • 数据持续到达的在线学习
  • 存储受限的嵌入式系统
  • 可以容忍一定精度损失的应用

7. 常见陷阱与错误 (Gotchas)

7.1 易被忽视的精度损失源

  1. 隐式的灾难性抵消 - 计算 $\mathbf{A}^T\mathbf{A}$ 形成法方程 - 两个接近矩阵的差：$\mathbf{A} - \mathbf{B}$ - 使用 $\log(1 + x)$ 而非 log1p(x)

  2. 维度诅咒引起的数值问题 - 高维空间中的距离计算 - 协方差矩阵在高维下趋向奇异 - 随机投影的概率保证vs实际精度

  3. 并行化引入的不一致性 - 不同的求和顺序 - 非原子操作的竞争条件 - 缓存一致性导致的微小差异

7.2 调试数值问题的系统方法

  1. 使用哨兵值： 在关键位置插入已知的测试向量 验证它们是否被正确处理

  2. 条件数监控： 定期估计和记录条件数 当条件数突变时触发警告

  3. 残差vs误差： 小残差 ≠ 小误差（病态问题） 大残差 ≠ 大误差（良态问题的舍入）

7.3 平台相关的陷阱

  1. 编译器优化的副作用 - -ffast-math 可能破坏IEEE合规性 - 循环向量化改变求和顺序 - 函数内联影响精度

  2. 硬件差异 - x86的80位扩展精度 - ARM的不同舍入模式 - GPU的非规范化数处理

  3. 库的不一致性 - 不同BLAS实现的数值差异 - MKL vs OpenBLAS vs参考实现 - 版本更新可能改变数值行为

8. 最佳实践检查清单

算法设计阶段

  • [ ] 条件数分析
  • 估计问题的固有条件数
  • 识别病态的子问题

  • 考虑重新表述以改善条件数

  • [ ] 算法选择

  • 优先选择后向稳定的算法
  • 考虑问题的特殊结构

  • 评估精度vs效率的权衡

  • [ ] 误差预算

  • 分配各步骤的误差容限
  • 识别误差放大的关键步骤
  • 设计误差补偿机制

实现阶段

  • [ ] 数值安全编程
  • 避免不必要的中间量计算（如 $\mathbf{A}^T\mathbf{A}$）
  • 使用稳定的基本操作（如 hypot, log1p）

  • 实现适当的缩放策略

  • [ ] 精度管理

  • 关键计算使用高精度
  • 实现混合精度策略

  • 保留误差估计信息

  • [ ] 鲁棒性增强

  • 添加数值健康检查
  • 实现优雅的降级策略
  • 提供诊断输出选项

测试阶段

  • [ ] 数值测试套件
  • 包含病态测试用例
  • 测试极端输入（很大/很小的数）

  • 验证误差界

  • [ ] 跨平台验证

  • 测试不同硬件架构
  • 验证不同编译选项

  • 检查并行化的数值一致性

  • [ ] 性能vs精度权衡

  • 基准测试不同精度级别
  • 评估算法的精度退化曲线
  • 确定最优工作点

部署阶段

  • [ ] 监控策略
  • 实时条件数估计
  • 异常检测机制

  • 性能退化预警

  • [ ] 文档要求

  • 明确数值假设和限制
  • 提供精度保证说明

  • 包含故障排除指南

  • [ ] 持续改进

  • 收集数值失败案例
  • 更新测试套件
  • 优化关键数值核心

特定场景补充

大规模计算

  • [ ] 考虑通信引起的误差
  • [ ] 实现检查点以应对数值失败
  • [ ] 使用分层算法减少误差累积

实时系统

  • [ ] 固定迭代次数避免不确定性
  • [ ] 预计算条件数阈值
  • [ ] 实现快速但略损精度的备选方案

机器学习应用

  • [ ] 梯度计算的数值稳定性
  • [ ] 防止权重更新中的数值崩溃
  • [ ] 批归一化等技术的数值考虑

延伸阅读与研究方向

  1. 自适应精度算法：根据问题条件数动态调整工作精度
  2. 概率数值方法：将舍入误差建模为随机变量
  3. 形式化验证：使用定理证明器验证数值算法
  4. 量子纠错码启发的数值方法：借鉴量子计算的错误处理
  5. 神经网络辅助的稳定性预测：使用ML预测数值失败模式

本附录提供了数值稳定性的实用指南。记住：在数值计算中，正确的算法只是起点，稳定的实现才是关键。当遇到数值问题时，系统地应用本附录的原则和技术，通常能找到解决方案。