JOVANA
Explore Library Glossary Getting Started Three Levels Fields How it works Mission
Join the mission
All guides

高斯消去法與 LU 分解

把學校裡那套消去變數的把戲,重新講成科學計算的主力演算法——再加上一個安靜的記帳洞見,把它變成一個可重複使用的分解 A = L U,讓你能對它一解再解。

長大成人的學校把戲

你早就會用手解一個小型方程組 A x = b:挑一個方程式,用它把某個變數從其他所有方程式裡消掉,重複到只剩一個變數,再倒著退回去求解。那套課本儀式有個成熟的名字,叫高斯消去法,而它是科學計算裡執行次數最多的單一演算法。本級的整個重點,是把它看成不是針對某一個 b 的一次性程序,而是一台你能乾淨地搖動、精確地計數、並且重複使用的機器——這正是 LU 分解登場的地方。

讓我們在一個小小的 3×3 方程組上搖動它,好讓所有零件都看得見。想像增廣陣列——矩陣 A 後面把右端 b 接上去當成最後一欄。演算法由左而右掃過去,一次處理一個主元欄。在每一欄它挑對角線上的元素當主元(pivot),然後對它下面的每一列,減去主元列的某個倍數,這個倍數恰好挑得讓主元正下方的元素變成零。那個乘數——你拿來縮放主元列的因子——正是你絕對不可以丟掉的那個數字。大多數人手算時,算出它就忘了它。記住它,就是 LU 的全部祕密。

[ 2  1  1 | 5 ]        [ 2  1  1 |  5 ]
[ 4  3  3 | 13]  -->   [ 0  1  1 |  3 ]   row2 -= 2*row1
[ 6  7  9 | 22]        [ 0  4  6 |  7 ]   row3 -= 3*row1

                       [ 2  1  1 |  5 ]
                 -->   [ 0  1  1 |  3 ]
                       [ 0  0  2 | -5 ]   row3 -= 4*row2

multipliers used:  2, 3, 4   (keep these!)
消去法的一次掃描:把每個主元下方都歸零。用到的乘數 2、3、4 就是 L 的元素。

兩個三角形免費掉出來

看看在例子裡消去法留下了什麼。變換後的矩陣是上三角的——對角線以下全是零——把它叫做 U。而我們被叮囑要留著的那些乘數(2、3、4),恰好填進一個對角線為 1 的下三角矩陣 L 的對角線以下位置。令人驚奇的事實是,這兩塊乘回去就得到原矩陣:A = L U。消去法不只解了一個方程組;它悄悄地把 A 分解成一個下三角矩陣乘上一個上三角矩陣。這兩個三角矩陣的乘積,就是 LU 分解

為什麼兩個三角形的乘積是這麼大的獎賞?因為三角形方程組解起來易如反掌。當 L 是下三角時,L y = b 可以從上往下直接讀出來——第一個方程式只有一個未知數,第二個有兩個但你已經知道第一個了,以此類推;這就是前向代入。當 U 是上三角時,U x = y 則用後向代入從下往上一層層剝開。所以一旦你有了 A = L U,解 A x = b 就拆成兩趟便宜的三角掃描。本級下一篇指南會專門細談這些代入到底怎麼運作,以及為什麼它們花費這麼少。

  1. 從 A x = b 出發並代入 A = L U,於是方程組讀作 (L U) x = b。把括號重組成 L (U x) = b。
  2. 把內層那塊命名為 y = U x。外層方程組現在是 L y = b,一個下三角方程組——用前向代入由上往下解它。
  3. 拿到 y 之後,解 U x = y,一個上三角方程組——用後向代入由下往上解它。你得到的 x 就滿足原來的 A x = b。

分解一次,求解多次

現在來看讓 LU 不只是高斯消去法漂亮改寫的那份回報。昂貴的部分——產生 L 與 U 的消去掃描——只碰 A,從不看 b。所以如果你必須對「同一個」矩陣 A 解許多個不同的右端 b_1、b_2、b_3、…,你只付一次消去成本,之後每個新的 b 都用一對極便宜的三角求解搞定。這就是分解一次、求解多次的想法,而它無所不在:一位結構工程師把上百種不同的載重情況推過同一個剛度矩陣;一段動畫釘住同一套物理,每一幀都重新求解。

數字本身說明了為什麼這很重要。對一個稠密的 n×n 矩陣,建出 L 與 U 的消去約需 (2/3) n^3 次運算——一件 O(n^3) 的工作,是主要的開銷。每次三角求解卻只約需 n^2 次運算,一件 O(n^2) 的工作。對 n = 1000 來說,那是約十億次運算去分解,對上約一百萬次去求解:第二個右端大約比第一個便宜一千倍。這個浮點運算計數不是學究氣;它是「跑一整夜」和「瞬間完成」之間的差別,也是為什麼沒有人會為每個 b 重跑消去。

天真的消去法在哪裡崩掉

上面的一切都假設了主元——我們拿來除的那個對角線元素——從來不是零。但你嘗試的第一個矩陣,對角線上就可能擺著一個零。這時乘數就成了除以零,演算法當場死亡。把那一列和下面某一列(在該位置有非零元素的)對調,就能立刻修好它,而這個列對調,正是選主元(pivoting)的種子。所以選主元不是後來才栓上去的可有可無的精修;有時候單純的消去法沒有它根本就走不下去。

更微妙的麻煩,是一個不是零、只是非常小的主元。除以一個小數字會讓乘數變得巨大,而巨大的乘數會放大浮點運算裡早已潛伏的捨入誤差——別忘了每個運算都會捨入,而且 0.1 沒有精確的二進位形式,所以機器從來不是在用你那些精確的數字。一個微小的主元可能把那些捨入誤差放大到計算出的答案變成胡言亂語,即使底下的問題本來規規矩矩。一個健全的問題,與一個把它毀掉的方法,之間的這道鴻溝,正是問題的條件數與演算法穩定性之間的差別。

解藥是部分選主元(partial pivoting):在某一欄消去之前,先掃描整欄找出絕對值最大的元素,把那一列換上來當主元。這讓每個乘數的大小最多是 1,從而把捨入控制住。有了部分選主元,高斯消去法在實務上就成了後向穩定的演算法——它回傳的是某個與你的問題只差一個捨入誤差的問題之精確答案,而那正是一個誠實的方法所能承諾的極限。下一篇指南為什麼我們要選主元會把穩定性的故事攤開;現在你只要記住這條規則:真正的求解器一律選主元,而它們實際算的是 P A = L U,其中 P 記錄了列對調。

結構,以及接下來的內容

那張 (2/3) n^3 的價格標籤,假設的是一個稠密矩陣,每個元素都可能非零。但從真實模型裡跑出來的矩陣往往是稀疏的——絕大多數是零,只有少數非零元素,分布成可預測的樣式。譬如格點上的熱方程,會給出一個非零元素緊貼對角線的帶狀矩陣。一個稀疏直接求解器從不儲存、也從不去碰那些零,所以一個稠密時要花 O(n^3) 的 LU 分解,對一條窄帶可以降到近乎 O(n)。一維的特例乾淨到有它自己的名字,叫 Thomas 演算法,是三對角方程組的一趟 O(n) 掃描。

不過有個誠實的陷阱。消去法可能把一個零變成非零——A 裡原本空著的位置,在 L 或 U 裡變成被佔據了。這叫填入(fill-in),而在一個排序糟糕的稀疏矩陣上,它可能把因子灌滿這麼多新的非零元素,以致稀疏的優勢蕩然無存。補救之道是先把行與列重新排序,挑一個讓填入盡量少的消去順序。這個填入與重排序的取捨,是稀疏直接求解的核心,而本級的最後一篇指南會認真回到成本與稀疏的主題。