第一卷悄悄做出的假设
在整个第一卷中,每一步高斯消元、每一个逆、每一个点积都是在精确算术下计算的。两数相加,得到精确的和。解一个方程组,得到精确的解。这个世界在计算机上并不存在。真实机器把数字存储在浮点运算中——一张由可表示值构成的有限网格,而每次运算都会把结果舍入到这张网格上。
最关键的数字是单位舍入误差 u,在 IEEE 双精度中约为 1.1e-16。模型很简单:当你计算 a op b 时,机器返回的是被一个不超过 u 的相对误差扰动后的真实结果,即 fl(a op b) = (a op b)(1 + d),其中 |d| <= u。单次运算几乎毫无损失。但十亿次运算串联起来,麻烦就可能累积。
# The classic shock, in any IEEE-double language: 0.1 + 0.2 == 0.3 -> False 0.1 + 0.2 -> 0.30000000000000004 # 0.1 is not representable exactly; it is stored as the # nearest grid point, and the tiny error survives the add. # Catastrophic cancellation: subtracting near-equal numbers a = 1.0000000001 b = 1.0000000000 a - b -> 1.00000008274e-10 (only ~1 good digit left) # The leading digits agreed and cancelled; what remains is # dominated by each operand's rounding error.
后向误差:评判答案的正确方式
问“我算出的 x 离真实 x 有多近?”(前向误差)往往是错误的问题,因为对于病态问题,即便完美的算法也无法接近。更好的问题是后向误差:我算出的 x 是否恰好是某个被轻微扰动的问题的精确解?若一个算法总是求解一个邻近的问题——仅被 O(u) 扰动——它就是后向稳定的。这是黄金标准,而且即便前向误差很大时也可达到。
这种区分令人释然。当真凶其实是矩阵中固有的 1e14 量级条件数时,你不再为算出离谱的 x 而责怪自己的代码。它还给你一个具体的检验:计算残差 r = b - A x_computed。若 ||r|| 相对于 ||A|| ||x|| 极小,那么你的求解器就是后向稳定的,仅此而已——哪怕 x 本身看起来与真相毫不相像。
廉价的保险:迭代精化
当后向稳定还不够、你还想把前向误差也压下来时,迭代精化常常几乎免费地实现这一点。其思路是:你已经分解了 A 来得到 x,于是复用这个分解去求解修正方程组的代价很低。用更高精度计算残差,求出修正量,加回去,再重复。
- 用稳定的分解先解一次 A x = b,得到 x_0。
- 构造残差 r = b - A x_0,最好用扩展精度,使它携带真实信息。
- 复用同一分解求解 A d = r——这很便宜,无需重新消元。
- 更新 x_1 = x_0 + d,然后循环直到残差不再缩小。