第 22 章 数值线性代数¶
前置：矩阵分解 (Ch10) · 范数与扰动 (Ch15) · 矩阵分析 (Ch14)
本章脉络：从理论解到数值解 \(\to\) 误差分析基础：绝对误差与相对误差 \(\to\) 浮点运算 (FLOPs) 与舍入误差 \(\to\) 后向误差分析 (Backward Error Analysis) \(\to\) 线性方程组的直接法（带有选主元的 LU、QR） \(\to\) 特征值算法：幂法、反幂法、QR 迭代 \(\to\) 稀疏矩阵技术 \(\to\) 迭代法基础：Jacobi, Gauss-Seidel, SOR \(\to\) 应用：超级计算机大规模仿真、工业级求解器 LAPACK 设计原理
延伸：数值线性代数是“能跑通的线性代数”；它研究的不是 \(Ax=b\) 的数学存在性，而是如何通过有限精度的计算机在 \(O(n^3)\) 或更低时间内得到最接近真相的解，是所有现代仿真软件的动力源泉
在理论数学中，我们假设实数具有无限精度。但在计算机中，一切数字都由有限位（如 64 位）表示。数值线性代数（Numerical Linear Algebra）专门研究在这种有限精度的约束下，如何设计高效且稳健的算法。一个优秀的数值算法不仅要快，更要在面对病态输入时保持结果的可靠性。本章将介绍这些构建现代科学计算大厦的底层算法与误差分析技术。
22.1 误差分析与浮点运算¶
定义 22.1 (相对误差)
若近似值为 \(\hat{x}\)，真实值为 \(x\)，则相对误差为 \(\frac{\|\hat{x} - x\|}{\|x\|}\)。 后向误差分析：与其研究结果偏离真实值多少，不如研究结果是哪一个“附近问题”的精确解。如果这个“附近问题”与原问题的距离极小，则算法是稳定的。
22.2 特征值计算：QR 迭代¶
定理 22.1 (特征值算法的局限)
Abel-Ruffini 定理证明不存在一般 5 次以上方程的求根公式。因此，所有的特征值算法（\(n \ge 5\)）本质上都是迭代法。 QR 迭代 是目前最通用的特征值算法： 1. \(A_k = Q_k R_k\) 2. \(A_{k+1} = R_k Q_k\) 随着迭代进行，\(A_k\) 逐渐收敛到上三角阵（特征值在对角线上）。
22.3 迭代求解线性方程组¶
核心迭代法
- Jacobi 迭代：基于 \(x_{k+1} = D^{-1}(b - (L+U)x_k)\)。
- Gauss-Seidel 迭代：利用已算出的最新分量，收敛通常比 Jacobi 快。
- 超松弛法 (SOR)：通过引入参数 \(\omega\) 进一步加速收敛。
练习题¶
1. [基础] 什么是 FLOPs？计算一个 \(n\) 维向量的点积需要多少 FLOPs？
参考答案
定义与计算： 1. FLOPs：Floating Point Operations，浮点运算次数。 2. 对于点积 \(\sum a_i b_i\)： - 需要 \(n\) 次乘法。 - 需要 \(n-1\) 次加法。 结论：总共约为 \(2n - 1\) 次运算。在大 \(O\) 记号下记为 \(O(n)\)。
2. [计算] 解释为什么 \(O(n^3)\) 的算法在 \(n=10^6\) 时不可接受。
参考答案
量级估算： 1. \(n = 10^6\) 时，\(n^3 = 10^{18}\) 次运算。 2. 假设一台每秒执行 \(10^9\) 次运算（1 GFLOPS）的计算机。 3. 需要时间 \(10^{18} / 10^9 = 10^9\) 秒。 4. \(10^9\) 秒约为 31.7 年。 结论：这就是为什么数值代数需要开发 \(O(n \log n)\) 或更高效算法的原因。
3. [稳定性] 在 LU 分解中，为什么要进行“选主元”（Pivoting）？
参考答案
理由： 如果主元 \(a_{ii}\) 接近 0，在计算乘子 \(l_{ji} = a_{ji}/a_{ii}\) 时会产生巨大的数值。 这会导致在后续减法中，原有的微小舍入误差被这个巨大乘子成倍放大，最终使 \(U\) 块的计算完全失真。选主元通过交换行，确保分母始终是当前列中绝对值最大的元素，从而将误差放大控制在最小范围。
4. [QR迭代] 若 \(A = Q R\)，为什么 \(R Q\) 的特征值与 \(A\) 相同？
参考答案
证明： 1. 已知 \(A = QR\)。 2. 由于 \(Q\) 是正交阵（酉阵），其逆矩阵为 \(Q^*\)。 3. 故 \(R = Q^* A\)。 4. \(RQ = (Q^* A) Q = Q^* A Q\)。 结论：\(RQ\) 与 \(A\) 是酉相似的。相似变换不改变矩阵的特征值。
5. [收敛性] Jacobi 迭代收敛的充要条件是什么？
参考答案
结论：谱半径 \(\rho(D^{-1}(L+U)) < 1\)。 在实践中，如果矩阵 \(A\) 是严格对角占优的，则 Jacobi 迭代保证收敛。
6. [稀疏矩阵] 对于 \(10^6 \times 10^6\) 的三对角矩阵，存储它需要多少内存？
参考答案
计算： 1. 三对角矩阵只有主对角线、上对角线和下对角线有值。 2. 非零元个数约为 \(3n = 3 \times 10^6\)。 3. 假设使用双精度浮点数（8 字节）：\(3 \times 10^6 \times 8 \approx 24\) MB。 对比：若存储为全矩阵，需要 \((10^6)^2 \times 8 = 8\) TB。 结论：稀疏存储技术是现代超大规模计算的生存基础。
7. [后向误差] 什么是算法的“后向稳定性”？
参考答案
如果算法在处理输入数据 \(x\) 时，得到的解 \(\hat{y}\) 恰好是另一个输入 \(\hat{x}\)（满足 \(\hat{x}\) 极其接近 \(x\)）的精确解，则称算法是后向稳定的。这意味着算法产生的误差并不比输入数据自带的噪声更严重。
8. [幂法] 简述幂法（Power Method）的作用及其局限。
参考答案
作用：用于求解矩阵的主特征值（模最大的特征值）。 局限： 1. 只能找主特征值。 2. 收敛速度取决于 \(\lambda_2 / \lambda_1\)。如果两个特征值很接近，收敛极慢。
9. [条件数] 若系统的条件数 \(\kappa(A) = 10^{12}\)，在 16 位精度的双精度运算下，解 \(x\) 可能有多少位是正确的？
参考答案
结论：约 4 位。 估算公式：有效位数 \(\approx\) 机器精度 - \(\log_{10} \kappa(A)\)。 计算：\(16 - 12 = 4\)。 这说明该问题高度敏感，计算结果只有前几位是可靠的。
10. [应用] 简述 LAPACK 和 BLAS 在数值计算中的关系。
参考答案
层次关系： 1. BLAS (Basic Linear Algebra Subprograms)：底层核心，负责最简单的向量加法、矩阵乘法。它针对特定硬件（如 CPU 缓存、指令集）进行了极致优化。 2. LAPACK (Linear Algebra Package)：上层架构，调用 BLAS 提供的原语来实现 LU、特征值、SVD 等复杂分解。 这种分层设计确保了数值软件在获得高层抽象的同时，拥有顶级的执行性能。
本章小结¶
数值线性代数是将数学真理转化为物理仿真的“最后一公里”：
- 稳定性的主宰：它揭示了算法的优劣不在于形式的优雅，而在于其对舍入误差的抑制能力，确立了后向误差分析的科学标准。
- 效率的极限：通过对 FLOPs 和内存布局的精细控制，数值代数突破了 \(O(n^3)\) 的屏障，使处理百万级维度的线性系统成为现实。
- 迭代的哲学：面对高次方程的不可求根性，迭代法成为了连接有限步骤与无限精度的桥梁，是描述连续动力学过程的最终手段。