災難性抵消與穩定的二次公式

陷阱：把兩個幾乎相等的數字相減

到現在你已經熟悉本級稍早的捨入運算模型：每個儲存的數字、每個基本運算都帶有大小至多為單位捨入 u 的微小相對誤差，在雙精度裡約為 10^(-16)。你或許會期待，只要每一步都精確到 16 位，答案就會維持 16 位精度。這個希望恰好會在一種常見情境下破滅，而且破滅得很戲劇化。它有個值得記住的名字：災難性抵消。

整個想法可以濃縮成一張小圖。假設兩個數字在前面幾位都一致——比如 a = 1.2345678 與 b = 1.2345600。它們各自都好好的。但把它們相減就得到 a - b = 0.0000078。七位一致的前導數字消失了；結果只保留了原本不同的兩位尾數。而那些尾數正是最不受捨入保護的位數——你先前看不見的雜訊，突然間變成了整個答案。

注意這裡的障眼法：減法本身算得幾乎完美。如果 a 與 b 正是你想要的精確數字，fl(a - b) 會精確到 u 之內，毫無問題。傷害早在 a 與 b 被儲存或被算出來時就造成了：它們各自帶有相對大小為 u 的誤差，換成絕對值大約是其量級（約 1.2）的 10^(-16) 倍。抵消之後答案只剩約 10^(-5)，於是同一個絕對誤差現在變成大約 10^(-16)/10^(-5) = 10^(-11) 的相對誤差——你十六位中的十一位沒了。減法並沒有製造誤差；它只是把原本看不見的捨入雜訊，從不顯眼提升為主宰。

條件 vs. 穩定性：誤差是誰的錯？

替兩種可能出錯的不同情況取名是有幫助的，因為災難性抵消是穩定性問題，而非條件問題。一個問題的條件數衡量真實答案對輸入微小變動的反應有多大——它屬於問題本身，跟你怎麼解它無關。一個演算法的數值穩定性衡量的，則是該演算法在這之上額外注入了多少誤差。稍早那條誠實的記帳規則：精度 ≈ 條件 × 穩定性。一個向後穩定的演算法，遇到條件惡劣的問題照樣會給出爛答案——而那裡沒有任何方法能救你。

災難性抵消正是「不穩定的演算法傷害一個條件良好的問題」的經典例子。從接近 1.23 的輸入算出 0.0000078，這件事本身是良性、條件良好的——真實 a 與 b 的微小相對變動，只會讓真實的差改變一個同樣微小的相對量。不穩定完全出在我們計算它的方式上：我們讓計算經過一個會放大輸入捨入誤差的減法。同樣的問題，不同的演算法，向前誤差卻天差地遠。

二次公式，以及它在哪裡崩潰

現在輪到那位著名的受害者。要解 a x^2 + b x + c = 0，你會自動伸手去拿每個學生都背過的公式：x = (-b ± sqrt(b^2 - 4 a c)) / (2 a)。在紙上它是精確的。在電腦上，它卻可能默默地交回一個幾乎沒有任何正確位數的根——而且失靈的情況一點也不稀奇，只要 b^2 遠大於 4 a c，使得 sqrt(b^2 - 4 a c) 非常接近 |b|，它就會發生。

看看那兩個分子。一個是 -b + sqrt(b^2 - 4 a c)，另一個是 -b - sqrt(b^2 - 4 a c)。當 sqrt(b^2 - 4 a c) ≈ |b| 時，這兩者中恰好有一個是兩個量級幾乎相等、符號相反的數的和——也就是一次幾乎相等者的相減——而另一個則是兩個同號數的安全和。安全的那個沒問題；另一個會抵消並弄丟大半位數。關鍵是，被毀掉的是量級較小的那個根，這很容易被忽略，因為較大的根看起來仍然完美。

Solve  x^2 + 200 x + 0.000015 = 0    (a=1, b=200, c=0.000015)

Exact roots:   x_1 = -199.99999992500...   x_2 = -7.5000000028e-8

Naive small root:  x = (-b + sqrt(b^2 - 4ac)) / (2a)
  sqrt(40000 - 0.00006) = 199.99999985
  -200 + 199.99999985  ->  -0.00000015   (two large near-equals subtract!)
  x ~ -7.5e-8   but only ~1-2 correct digits survive

Stable small root (use the good root + product of roots):
  x_large = (-b - sqrt(b^2 - 4ac)) / (2a)   (safe: same-sign sum)
  x_small = c / (a * x_large)               (x_1 * x_2 = c/a, no subtraction)
  x_small ~ -7.5000000028e-8   (full accuracy)

天真版的小根去相減兩個幾乎相等的數，弄丟了大半位數；改走 x_1 x_2 = c/a 則完全避開了那個減法。

那一行修正，以及它為何有效

解方是一個小小的代數重排，它從不去相減幾乎相等者。先算出那個分子為安全同號和的根——選擇符號，使你把 -b 與平方根以一致的符號相加，也就是 q = -(b + sign(b)·sqrt(b^2 - 4 a c)) / 2。這可靠地給出一個根 x_1 = q / a。接著，那個危險的根不靠會抵消的公式來求，而是用根的乘積：對 a x^2 + b x + c = 0 而言，x_1 · x_2 = c / a。於是 x_2 = c / q。一個除法，沒有減法——抵消乾脆就不會發生。

計算判別式 d = b^2 - 4 a c。（若 d < 0 則根為複數；若 d 極小，兩根幾乎重合，那是問題本身另一種真正的條件惡劣，不是這個改寫能修好的。）
構造 q = -(b + sign(b)·sqrt(d)) / 2。透過讓符號與 b 一致，你相加的是兩個同號量——這是不會抵消的安全運算。
把兩個根回傳為 x_1 = q / a 與 x_2 = c / q。兩者都只用乘法與除法算出；沒有任何一條路徑去相減兩個幾乎相等的數。

這就是數值穩定性在實務中的樣子：同樣的數學問題、同樣的輸入，但走一條不會放大那不可避免之輸入捨入的算術路徑。天真版與穩定版解的是同一個、條件良好的問題；只有穩定版守住了完整的約 16 位，因為它的向後誤差很小——算出的根，恰好是某個係數與你的僅相差一根毛髮的二次式的精確根。這正是我們追求的務實目標：不是零誤差，而是某個鄰近問題的精確答案。

同一招用在其他所有地方

一旦你練出了眼力，就會在數值工作裡到處看見抵消的陷阱，而同一個「改寫代數式」的反射動作就能化解它們。幾個你在這條階梯上會再遇到的：對小的 x 計算 1 - cos(x)（改寫成 2 sin(x/2)^2，或用專門的函式）；變異數公式 E[x^2] - E[x]^2（改用兩遍法或 Welford 形式）；對小的 x 計算 exp(x) - 1（用 expm1）；以及用極小的 h 算差商 (f(x+h) - f(x)) / h，其分子去相減兩個幾乎相等的值，結果撞上捨入下限——這正是為什麼把 h 縮小並不會讓有限差分導數持續變好。

帶兩個誠實的提醒繼續前進。第一，改寫治的是不穩定，從不治條件惡劣：若底層問題本身真的敏感——一個有真正重根的二次式、一個條件數接近 10^8 的矩陣——那麼再聰明的代數也救不回問題自己都釘不住的位數；預期會在 16 位中損失大約 8 位。第二，每個浮點答案仍然是受機器 epsilon支配的近似值；一個穩定演算法的目標不是精確，而是在問題已經強加於你的誤差之上，盡可能少加一點。