線性代數/主題：計算精度

高斯消元法非常適合計算機化。下面的程式碼說明了這一點。它在一個$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.