為什麼三角形是最容易的情形
在上一篇,你看著高斯消去法把一個稠密矩陣 A 輾成一個分解 A = L U,其中 L 是下三角(對角線以上全為零),U 是上三角(對角線以下全為零)。那看起來工程浩大,事實也是如此——大約是 n^3 的三分之二次運算。很自然的問題是:何苦呢?一個三角系統憑什麼就比原本那團亂麻更好?這篇就是答案,而且是個令人滿意的答案:三角系統幾乎可以免費求解,在一趟整齊的掃描裡完成,而這趟掃描,才是真正把你的 x 交到手上的東西。
畫面是這樣的。一般的系統 A x = b 把每個未知數都跟其他所有未知數綁在一起:第 1 列同時提到 x_1 到 x_n,你不先解開其餘的就無法釘住任何單一未知數。三角系統解開了這個結。在一個上三角系統 U x = b 裡,最後一條方程只含一個未知數——它寫成「(某數)乘以 x_n =(某數)」,所以 x_n 立刻靠一次除法掉出來。一旦知道了 x_n,倒數第二條方程就只剩 x_{n-1} 還未知,於是它也能被除出來。你一階一階地爬這道樓梯,每一階都恰好有一個新的未知數獨自等著你。
回代法,一步步走過
那個從底部往上爬樓梯的演算法叫做回代法(back substitution),它求解上三角的那一半 U x = b。這名字是字面的意思:你先找到最後一個未知數,然後把它代回上面那條方程,那裡就空出一個新的未知數,再重複。我們親手做一個小小的 3×3,好讓這個動作清清楚楚。
Solve U x = b :
2 x1 + 1 x2 - 1 x3 = 3
3 x2 + 2 x3 = 5
4 x3 = 8
Row 3: 4 x3 = 8 -> x3 = 2
Row 2: 3 x2 + 2(2) = 5 -> 3 x2 = 1 -> x2 = 1/3
Row 1: 2 x1 + (1/3) - 2 = 3 -> 2 x1 = 3 - 1/3 + 2 = 14/3 -> x1 = 7/3- 從最底列開始,它只含一個未知數 x_n;把它的右端項除以對角線元素,就得到 x_n。(那個對角線元素就是主元——它不能是零,而這正是下一篇講樞軸選取的篇章所要保證的事。)
- 往上移一列。它下面的每個未知數現在都已知,所以把那些值代進去;剩下的又是「一條方程、一個未知數」。
- 把已知的那些項從右端項裡減掉,再除以這一列的對角線元素,就放出下一個未知數。
- 重複到你抵達最頂列,最後一個未知數 x_1 掉出來為止。你現在從上到下擁有了整個向量 x。
前代法是它的鏡像。它求解一個下三角系統 L y = b,只是現在那條容易的方程是第一條(它只提到 y_1),所以你從頂端開始、往下走。動作一模一樣,只是上下顛倒。我們之所以對兩者都在意,是因為 A = L U 把一次困難的求解,拆成兩次容易的求解,下一節會把這點講精確。
兩次三角求解,取代一次困難的求解
現在我們可以看清楚,為什麼那番分解的汗水沒有白流。一旦有了 A = L U,要解 A x = b,你就再也不碰 A。你做代換:A x = b 變成 L U x = b。把中間那塊命名為 y = U x。於是 L y = b 是下三角的——用前代法解它得到 y。接著 U x = y 是上三角的——用回代法解它得到 x。一個真正困難的問題,就這樣依序變成了兩個真正容易的問題。
A x = b with A = L U step 1 (forward substitution): L y = b -> y [lower triangular] step 2 (back substitution): U x = y -> x [upper triangular]
把成本數一數,策略就不證自明了。一次三角求解大致把三角形裡每個元素碰一次:對一個 n×n 的三角形,那大約是 n^2/2 次乘加運算,所以兩次合起來約 n^2——一個 O(n^2) 的運算量。拿它跟「建出 L 與 U」那場消去的 O(n^3) 相比。當 n = 1000,n^3 比 n^2 大上一千倍。昂貴的步驟是分解;而真正產出 x 的那兩趟三角求解,相比之下便宜得很。
這道成本鴻溝,正是「把 LU 分解存下來、留著用」的全部理由。假設你必須對同一個 A、卻有許多不同的右端項 b 求解 A x = b——同一座橋上不同的載重、同一個電路在不同的日子。你把 A 分解成 L U 恰好一次,那筆 O(n^3) 的帳只付一遍,之後每個新的 b 只花一對 O(n^2) 的三角掃描。這就是著名的分解一次、求解多次模式,而三角求解正是讓它划算的、便宜的後半段。
這趟掃描可信嗎?穩定性與主元
一個合理的擔憂:我們在浮點數裡做一長串的除法與減法,而那裡的算術並不精確——0.1 沒有精確的二進位形式,加法甚至不滿足結合律。當我們往上爬時,誤差會不會危險地堆積起來?令人安心而誠實的答案是:不會,至少對三角求解而言。回代法與前代法是向後穩定的:算出來的 x,恰好是某個三角系統的精確解,而那個系統的元素跟真正的元素只差一根毫毛,數量級就在單位捨入誤差上下。這正是數值計算裡務實的黃金標準——不是一個剛好正確的答案,而是「一個離你的問題只有一髮之遙的問題」的剛好正確的答案。
這趟掃描仰賴一個絕對的要求:你拿來當除數的每個對角線元素——每個主元(pivot)——都必須非零。如果 U 的對角線上有一個零,那麼矩陣是奇異的,根本不存在唯一的 x 可找;那次除法會炸成無窮大或 NaN,而那正是演算法在誠實地回報:這個問題沒有答案。即使主元只是很小(不是剛好為零)也是個警訊,因為除以一個接近零的數會放大誤差。讓主元安全地遠離零,恰恰是部分樞軸選取的工作,也就是緊接著下一篇的主題——好讓你餵進回代法的那些三角因子,一開始就有性質良好的對角線。
為什麼這也正是你永不求逆的理由
初學者常在程式裡寫 x = A^{-1} b,照著課本公式逐字照抄。三角求解正是你不該這麼做的理由。要算出 A^{-1},你得求解 A X = I——那是 n 個各自獨立的右端項(單位矩陣的各行),所以是在分解之上再加 n 對三角求解,比你需要的工作量更多。然後再把 A^{-1} 乘上 b,又是另一個 O(n^2) 的步驟,添上它自己的捨入誤差。直接用一次前代加一次回代來求解,既更便宜又更準確,因為你從頭到尾都沒有建出、也沒存下那個易出錯的反矩陣。
最後一個朝前看的連結。正因為三角求解既便宜又穩定,它成了遠超出本級的、可重複使用的積木。當矩陣是稀疏或帶狀的——比方一個三對角系統,其中的托馬斯演算法不過就是前代再回代、被精簡到只剩三條對角線——同樣這趟掃描就以 O(n) 而非 O(n^2) 跑完。而在後面的篇章裡,前置條件器與迭代精化會一再倚賴一次快速的三角求解作為它們的內層步驟。你剛學會的這趟不起眼的爬樓梯,是整個科學計算裡最常被用到的核心運算之一。