线性代数/主题：计算精度

高斯消元法非常适合计算机化。下面的代码说明了这一点。它在一个$n\!\times \!n$ 矩阵上操作a, 首先用第一行进行主元运算，然后用第二行进行主元运算，等等。

for (pivot_row = 1; pivot_row <= n - 1; pivot_row++) {
    for (row_below = pivot_row + 1; row_below <= n; row_below++) {
        multiplier = a[row_below, pivot_row] / a[pivot_row, pivot_row];
        for (col = pivot_row; col <= n; col++) {
            a[row_below, col] -= multiplier * a[pivot_row, col];
        }
    }
}

(此代码是用 C 语言编写的。以下是简要的翻译。循环结构for (pivot_row = 1; pivot_row <= n - 1; pivot_row++) { ... }设置pivot_row为 1，然后在pivot_row小于或等于 $n-1$ 的情况下进行迭代，每次迭代都用 “pivot_row” 操作将++增加 1。另一个不明显的结构是 “-=” 在最内层的循环中相当于a[row_below, col] =- multiplier * a[pivot_row, col] + a[row_below, col]}操作。）

虽然这段代码提供了对高斯消元法如何机械化的快速了解，但它并不准备使用。它在很多方面都很天真。最明显的是，它假设在pivot_row, pivot_row位置上始终可以找到一个非零数用作主元元素。为了使其实用，重写这段代码的一种方法是涵盖在该位置找到零导致行交换的情况，或者得出矩阵是奇异的结论。

添加一些if$\cdots$ 语句来涵盖这些情况并不难，但我们将考虑这段代码在哪些方面天真。由于计算机依赖有限精度浮点运算，因此存在陷阱。

例如，我们已经看到，我们必须将奇异系统作为单独的情况处理。但是，几乎奇异的系统也需要小心。考虑以下系统。

${\begin{array}{*{2}{rc}r}x&+&2y&=&3\\1.000\,000\,01x&+&2y&=&3.000\,000\,01\end{array}}$

通过观察，我们得到解$x=1$ 且 $y=1$。但计算机很难做到。一台将实数表示为八位有效数字（通常称为单精度）的计算机将内部表示第二个方程为$1.000\,000\,0x+2y=3.000\,000\,0$，丢失了第九位的数字。这台计算机不会报告正确的解，而是报告一个与之相差甚远的解——这台计算机认为该系统是奇异的，因为这两个方程在内部被表示为相等。

为了直观地了解计算机如何得到一个相差甚远的解，我们可以绘制系统的图形。

在这个图形的尺度上，这两条线无法区分开来。该系统几乎是奇异的，因为这两条线几乎是同一条线。接近奇异性赋予该系统一个性质，即系统的小变化会导致解的大变化；例如，将$3.000\,000\,01$ 更改为$3.000\,000\,03$ 将交叉点从 $(1,1)$ 更改为 $(3,0)$。该系统根据第九位数字发生了根本性的变化，这解释了为什么八位计算机难以处理。对输入值中的不准确或不确定性非常敏感的问题称为病态。

以上示例提供了一种系统在计算机上难以求解的方式。它具有优势，即近似等线的图像可以让人们记住数值困难出现的一种方式。不幸的是，当我们想要解决一些大型系统时，这种洞察力并不是很有用。我们通常无法希望理解任意大型系统的几何形状。此外，除了某些线性曲面之间的角度很小之外，计算机的结果还可能以其他方式不可靠。

例如，考虑下面的系统，来自（Hamming 1971）。

${\begin{array}{*{2}{rc}r}0.001x&+&y&=&1\\x&-&y&=&0\end{array}}\qquad (*)$

第二个等式给出 $x=y$，因此 $x=y=1/1.001$，因此两个变量的值都略小于 $1$。使用两位数的计算机以这种方式在内部表示该系统（我们将用两位数浮点算术来做这个例子，但用八位数做一个类似的例子很容易）。

${\begin{array}{*{2}{rc}r}(1.0\times 10^{-3})x&+&(1.0\times 10^{0})y&=&1.0\times 10^{0}\\(1.0\times 10^{0})x&-&(1.0\times 10^{0})y&=&0.0\times 10^{0}\end{array}}$

计算机的行化简步骤 $-1000\rho _{1}+\rho _{2}$ 会产生一个第二个等式 $-1001y=-999$，计算机将其四舍五入到两位数，得到 $(-1.0\times 10^{3})y=-1.0\times 10^{3}$。然后计算机从第二个等式中确定 $y=1$，并从第一个等式中确定 $x=0$。这个 $y$ 值相当好，但 $x$ 很糟糕。因此，另一个导致输出不可靠的原因是浮点算术与对小枢轴的依赖的混合。

有经验的程序员可能会说，我们应该使用双精度，其中保留十六个有效数字。这确实可以解决很多问题。但是，作为一种普遍的方法，它也有一些困难。一方面，双精度比单精度慢（在 '486 芯片上，单精度乘法需要 11 个时钟周期，而双精度乘法需要 14 个时钟周期（Microsoft 1993）），并且内存需求是单精度的两倍。因此，尝试以双精度进行所有计算并不实际。 此外，上面的系统显然可以调整，以便在第十七位出现相同的麻烦，因此双精度并不能解决所有问题。我们需要一个策略来最大限度地减少在计算机上求解系统时出现的数值问题，以及一些关于如何确定报告的解决方案的可信度的指导。

数学家已经对如何获得最可靠的结果进行了深入研究。对上述朴素代码的一个基本改进是不简单地取pivot_row中pivot_row位置的元素作为主元，而是查看pivot_row列中pivot_row行以下的所有元素，并取最有可能给出可靠结果的元素（例如，取一个不太小的元素）。这种策略称为部分主元法。例如，要解上述有问题的系统($*$)，我们首先查看两个方程，寻找最佳第一个主元，并取第二个方程中的$1$，因为它更有可能给出良好的结果。然后，$-.001\rho _{2}+\rho _{1}$的主元步骤得到第一个方程为$1.001y=1$，计算机将表示为$(1.0\times 10^{0})y=1.0\times 10^{0}$，得出结论$y=1$，然后回代得到$x=1$，这两个结果都接近正确值。上述代码可以适应此目的。

for (pivot_row = 1; pivot_row <= n - 1; pivot_row++) {
    /* Find the largest pivot in this column (in row max). */
    max = pivot_row;
    for (row_below = pivot_row + 1; pivot_row <= n; row_below++) {
        if (abs(a[row_below, pivot_row]) > abs(a[max, row_below]))
            max = row_below;
     }

    /* Swap rows to move that pivot entry up. */
    for (col = pivot_row; col <= n; col++) {
        temp = a[pivot_row, col];
        a[pivot_row, col] = a[max, col];
        a[max, col] = temp;
     }

     /* Proceed as before. */
     for (row_below = pivot_row + 1; row_below <= n; row_below++) {
         multiplier = a[row_below, pivot_row] / a[pivot_row, pivot_row];
         for (col = pivot_row; col <= n; col++) {
             a[row_below, col] -= multiplier * a[pivot_row, col];
         }
     }
}

对实现高斯消元法的最佳方法的全面分析超出了本书的范围（见(Wilkinson 1965)），但大多数专家推荐的方法是上述代码的一种变体，它首先在候选主元中找到最佳主元，然后将其缩放到一个不太可能出现问题的数字。这称为带缩放的部分主元法。

除了返回一个可能可靠的结果外，大多数完成良好的代码还会返回一个称为条件数的数字，它描述了输入数字的不确定性被放大成为返回结果中的不精确度的倍数（见(Rice 1993)）。

本次讨论的教训是，仅仅因为高斯消元法在理论上总是有效的，仅仅因为计算机代码正确地实现了该方法，仅仅因为答案出现在绿条纸上，并不意味着答案是可靠的。在实践中，总是使用一个专家们已经努力克服可能出现问题的软件包。

练习

问题 1

使用两位小数，将$253$和$2/3$相加。

问题 2

这个求两条直线交点的問題与上述例子形成对比。

                 $\displaystyle {\begin{array}{*{2}{rc}r}x&+&2y&=&3\\3x&-&2y&=&1\end{array}}$

说明在这个系统中，数字的微小变化只会导致解的微小变化，方法是将底部方程中的常数改为$1.008$，并求解。将其与未改变的系统的解进行比较。

问题 3

用手解这个系统(Rice 1993)。

${\begin{array}{*{2}{rc}r}0.000\,3x&+&1.556y&=&1.569\\0.345\,4x&-&2.346y&=&1.018\end{array}}$

1. 用手精确求解。
2. 在每一步都四舍五入到四位有效数字，求解。

问题 4

计算机内部的舍入通常会影响结果。假设你的机器有八位有效数字。

1. 证明机器将计算 $(2/3)+((2/3)-(1/3))$ 与 $((2/3)+(2/3))-(1/3)$ 不相等。因此，计算机算术不具有结合性。
2. 比较计算机版本中的 $(1/3)x+y=0$ 与 $(2/3)x+2y=0$。第一个方程的两倍与第二个方程是否相同？

问题 5

病态性不仅取决于系数矩阵。这个例子 (Hamming 1971) 表明它可能源于系统左右两侧的交互作用。令 $\varepsilon$ 是一个小的实数。

${\begin{array}{*{3}{rc}r}3x&+&2y&+&z&=&6\\2x&+&2\varepsilon y&+&2\varepsilon z&=&2+4\varepsilon \\x&+&2\varepsilon y&-&\varepsilon z&=&1+\varepsilon \end{array}}$

1. 用手解该系统。注意，$\varepsilon$ 能够消去仅仅因为右侧和左侧的整数部分都完全抵消了。
2. 用手解该系统，四舍五入到小数点后两位，并且 $\varepsilon =0.001$。

解答

参考资料

• Hamming, Richard W. (1971), Introduction to Applied Numerical Analysis, Hemisphere Publishing.
• Rice, John R. (1993), Numerical Methods, Software, and Analysis, Academic Press.
• Wilkinson, J. H. (1965), The Algebraic Eigenvalue Problem, Oxford University Press.
• Microsoft (1993), Microsoft Programmers Reference, Microsoft Press.