数值计算的稳定性与有效数字位数

投稿时间：2013-09-18 22:39

牛津大学计算中心主任、英国皇家学会院士、美国工业与应用数学会（SIAM）主席Nick Trfethen 教授2012年11月的《SIAM News》上撰文链接内容以短小篇幅简明地阐述了数值线性代数方向两个未解之谜：1. 如何快速求逆矩阵 2. 高斯消元法在实际应用中的稳定性。

第一个问题的原文如下：

The first problem is, can an n × n linear system of equations Ax = b be solved in O(n^2+ε) operations, where ε is arbitrarily small? In other words, could there be a “fast matrix inverse,” like the fast Fourier transform?

第二个问题的原文如下：

The second big unsolved problem is also in the area of linear algebra. Why is Gaussian elimination numerically stable in practice, even though it’s unstable in the worst case? In other words, why does Gaussian elimination work?

数值计算稳定性的表现是：某特定问题的计算，在不同的数字计算机上的结果是不同的。则这个算法是不稳定的。稳定的基本含义是，某方法得出的结果和其真实值很接近（最好一样），且受不同具体计算机的影响很小。

有效数字位数有限造成的舍入误差（round-off error）是数字计算机稳定性的本质根源。截断误差、数据误差，是数字计算机以外的原因。

在代数运算里，有的计算引起的有效数字位数变化小，有点则大。“引起有效数字位数变化大的计算”是形成数值稳定性的主要原因。例如，常见的加法、减法的结果的有效数字位数变化慢。两个同长度的加法，引起有效数字位数的增加速度为O(lgn)，这里n为加法的次数。乘法、除法则引起O(n)的有效数字位数变化。

其中除法最值得关注：小绝对值的分母（特别是无理数分母）导致较大（亦即较长）的除法计算结果。当有效数字位数有限时，造成舍入误差会被明显放大。小绝对值的分母，涉及矩阵时表现为矩阵的奇异性。

降低或限制乘法、除法（包括矩阵形式的）引起的有效数字位数变化，是增加稳定性的核心路径之一。

数值形式的高斯消元法，不直接涉及矩阵计算，对N元方程，只涉及O(N³)的通常四则运算。高斯消元法里的乘法，几乎都是未知量和“能表示的常数”之间进行；不涉及两个或多个未知量之间的乘除法（即不是用3次以上的连续乘除法）。这样，对于同长度的系数和未知量，乘法除法造成的有效数字位数变化是O(1)的，即固定位数的。这样，整个高斯消元法造成的有效数字位数的变化速度为O(lgN)，即其中的加减法引起。

类似地，按照定义进行矩阵乘法，也是一种比较稳定的数值计算方法。对于元素同长度的N阶数值方阵，其中的乘除法引起的有效数字位数变化是O(1)的；加减法引起O(lgN)的有效数字位数变化。以Strassen算法为代表的O(N^2+ε)类矩阵乘法（这里实数0≤ε<1），核心是“把多个乘法合并在一起”计算（例如先将两个矩阵元素相加，再做乘法Q1=(a11+a22)*(b11+b22)：这样直接增大了乘法必须使用的有效数字位数。在特定有效数字位数的计算机上，会造成更大的截断误差）。一般地，O(N^2+ε)类矩阵乘法中的ε越小，其中使用的乘法必须的有效数字位数会越大，从而数值计算的稳定性会越差。

这类方法的极限是1990年陈道琦、谢友才、应文隆在科学通报发表的《关于矩阵乘法的一个最佳算法》。该方法只用一次乘法，但需要Ω(n³×lgn)的有效数字位数。

当参加加减乘除运算的元素长度不同时，一般会进一步增大有效数字位数的变化，即产生更不稳定的计算结果。

参考文献：

[1] Nick Trefethen. The Smart Money Is on Numerical Analysts [J]. SIAM News, Volume 45, Number 9, November 2012.

[2] 美国SIAM主席、皇家学会院士提的两个数值代数问题（1）[EB/OL]. http://www.mysanco.com/wenda/index.php?class=discuss&action=question_item&questionid=1097

[3] 杨正瓴，王渭巍，曹东波，张军，陈曦. 短期负荷预测的Ensemble混沌预测方法[J]. 电力系统自动化, 2007, 31(23): 34-37.

[4] Volker Strassen. Gaussian Elimination is Not Optimal [J]. Numerische Mathematik, 1969, 13(4): 354-356.

[5] 陈道琦，谢友才，应文隆. 关于矩阵乘法的一个最佳算法[J]. 科学通报, 1990, 35(3): 161-161.