淺談數值線性代數
賴玉玲
一 . 介紹
什麼是數值線性代數呢? 從線性代數 這一學科中, 我們學到了琳瑯滿目的方法來 解決線性代數問題, 這些方法都有一個共同 的特性, 計算量隨著矩陣維數 (dimension) 的增大而可能增加到非人力可以勝任, 以至 於必須訴諸電子計算機的幫助。 由於電腦的 記憶體, 運算速度及算術誤差等因素, 我們 不能單純的利用電腦來執行我們在線性代數 學科中所學的方法, 因此我們得針對電腦的 特性, 並藉助發展了好幾世紀的線性代數理 論, 設計出快速有效及準確的演算法 (algo- rithms), 這也就是數值線性代數這一學科的 主要研究課題。
眾所周知, 數學算子中最簡單的就是線 性算子, 它也是在應用數學領域中最常見的 模型, 所以其衍生出的線性代數問題值得探 討。 大體來說, 數值線性代數比較常處理其中 的兩大基本問題
(1) 找線性系統 Ax = b 的解
(2) 特徵值問題 (找特徵值 λ 及特徵向 量 x 使得 Ax = λx)
根據估計約有 75% 的應用數學問題會衍生 出以上的第一個問題, 其中又以偏微分方程
的數值解為最常見, 其次為函數的線性插值 (linear interpolation), 線性逼近 (lin- ear approximation) 及最小方差 (least square) 等問題, 所以在此我們將對第一個 問題做較深入的探討, 希望藉此可引導讀者 們進入數值線性代數領域。
求線性系統解的數值方法眾多, 大致可 粗略的分成兩大類 (1) 直接法 (2) 疊代法, 我們將於下列兩節分別詳述之。
二 . 直接法
直接法的特性在於如果不考慮電腦計算 的誤差, 則在有限步驟內就可得到正確解, 其 中最常被採用的即是大家所熟知的高斯消去 法 (Gaussian elimination) 及其變形 (vari- ants), 其主要原因是它的觀念淺顯易懂及它 的演算法可輕易地翻譯成電腦程式。
為了詳細描述它的觀念, 我們考慮下列 的線性系統
a11x1+ a12x2+ · · · + a1nxn= b1
a21x1+ a22x2+ · · · + a2nxn= b2
... ...
an1x1+ an2x2+ · · · + annxn= bn
(1) 假如 a11 6= 0, 則第 (i) 式減去第 (1) 式乘以
1
ai1
a11, 我們得到下列維數 (n − 1) 的線性系統
a(2)22x2+ a(2)23x3+ · · · + a(2)2nxn = b(2)2
a(2)32x2+ a(2)33x3+ · · · + a(2)3nxn = b(2)3
... ...
a(2)n2x2+ a(2)n3x3+ · · · + a(2)nnxn = b(2)n (2) 其中 a(2)ij = aij−aai1
11a1j, b(2)i = bi−aai1
11b1, i, j = 2, . . . n。 很顯然地, 變數 x1 被消去 了, 整個問題被簡化成 (n − 1) 的線性系統, 此時如果 a(2)22 6= 0, 則可重覆上述的過程, 更 甚者, 如果 a(k)kk 6= 0, k = 2, 3, . . . (n − 1), 則可依序消法變數 x2, x3, . . . xn−1直到只剩 下一個式子及一個變數。 如果把每一步驟中 的第一個式子收集起來, 我們得到下列三角 (triangular) 系統
a(1)11x1+ a(1)12x2+ · · · + a(1)1nxn = b(1)1 a(2)22x2+ · · · + a(2)2nxn = b(2)2
. .. ...
a(n)nnxn = b(n)n (3) 其中為了符號的一致性, 以 a(1)1j 及 b(1)1 分別 代表 aij 及 b1。 以上的過程就是所謂的高斯 消去法, 我們以下列的演算法做一個總結
for k = 1 to n − 1 for i = k + 1 to n
c := aik
akk
bi := bi− c bk
for j = k + 1 to n aij := aij − c akj
end
end end
利 用高 斯 消 去 法, 原 來 的線 性 系 統 Ax = b 被簡化成三角線性系統 Ux = ˜b, 其中 U 是一上三角矩陣。 新的問題 Ux = ˜b 十分易解, 從方程組 (3) 可輕易觀察出 xn = b(n)n /a(n)nn, 一旦 xn 值知道了, 代入 第 (n − 1) 式就可得到 xn−1 的值, 依此類 推, 可陸續得到 xn−2, xn−3, . . . x1 的值, 此 法被稱為回置代入法 (back substitution)。
高斯消去法雖是一淺顯易懂的方法, 但 是從它的推演過程中, 我們不難發現它並不 是一個進行無礙的方法, 因為途中若某個值 a(k)kk = 0 時 ( a(k)kk 被稱為 pivot), 則高斯 消去法就被迫中斷。 利用線性代數的理論, 不 考慮算術誤差的情況下, 可推導出某些矩陣 可以高枕無憂的使用高斯消去法, 然而由於 無法避免電腦計算的誤差, 所以對那些理論 上可使用高斯消去法的矩陣還是有可能遇到 pivot 為零的情況而被迫中斷高斯消去法的 使用, 因此我們不詳述那些矩陣理論上可無 礙的使用高斯消去法, 然而我們將看是否有 補救之道呢? 首先我們必須注意到高斯消去 法中斷並不表示此線性系統無解, 例如考慮
x1+ x2+ x3 = 6 x1+ x2+ 2x3 = 9 x1− x2+ 4x3 = 11
(4)
其解為 x1 = 1, x2 = 2 及 x3 = 3。 應用一 步高斯消去法後, 得到
(
0x2+ x3 = 3−2x2+ 3x3 = 5 (5) 很顯然地, 因為 a(2)22 = 0, 無法再使用高斯消 去法了。 然而我們可以藉著交換行、 交換列或
交換行及列來避免 pivot 為零的情況, 這種 方式通常被稱為 pivoting。 在更深入探討此 法前, 首先我們要提醒讀者注意以下兩項事 實: (1) 兩列互換的同時, 右邊向量 b 的兩 個相對應元素也要互換, (2) 兩行互換時不要 忘了相對應的兩個變數也要互換。
假如我們處理的矩陣是可逆的, 理論上 高斯消去法配合 pivoting 的技巧絕對可以 進行無礙地把解找到, 但實際上, 由於計算的 誤差, 我們無可避免地去考慮如何 pivoting 才能真正如理論上所說順暢地把解求出, 更 甚者, 如何 pivoting 才能增加解的準確性, 有許多不同的 pivoting 技巧可以達到上述 的目的, 在此我們列舉兩個最常被使用的
1. Complete pivoting: 在所處理矩陣的 所有元素中, 找一個絕對值最大的元素, 再利用行、 列交換把它搬到 pivot 的位 置。 以數學語言來說, 如果 r 及 s 使得
|a(k)rs | = maxi,j=k...n|a(k)ij |, 則第 r 列與 s 行分別與第 k 列及行互換。
2. Partial pivoting: 在所處理矩陣中的 第一行元素中, 找一個絕對值最大的元 素, 利用列交換把它搬到 pivot 的位 置, 換言之, 如果 r 使得 |a(k)rk| = maxi=k...n|a(k)ik |, 則第 r 列與第 k 列互 換。
因在 partial pivoting 的技巧中, 列交 換比較常被使用, 所以我們對其略做了些描 述, 實際上, 上述的觀念也可以考慮行交換來 代替之。 為了使讀者更易了解 pivoting 的 技巧, 我們考慮前述的例子, 利用 complete
pivoting, 我們得到
(
3x3− 2x2 = 5x3+ 0x2 = −1 (6) 雖然線性系統不變, 但是處理的矩陣 A2, 變 成了
"
3 −2 1 0
#
以 3 為 pivot, 相對地, 如 果考慮 partial pivoting, 則矩陣 A2 變成 了
"
−2 3 0 1
#
以 −2 為 pivot。
理論上高斯消去法配合 complete piv- oting 所得到的解較配合 partial pivoting 而得的解來得精確, 但在眾多的例子中, 數值 結果顯示配合 partial pivoting 的技巧所得 到的解跟實際的解差異並不大, 而且從 com- plete pivoting 及 partial pivoting 的定義 中, 可輕易看出 partial pivoting 所需的工 作量遠少於 complete pivoting 的工作量, 所以一般的情況下, 除非想得到非常精確的 解, partial pivoting 技巧還是較常被建議 的。
三 . 疊代法
大型電路 (large circuits) 分析、 氣象 學、 地震學, 及太空科技等推演出的偏微分方 程, 經離散化後所導出的線性系統, 一般來說 都有成千上萬的未知數 (亦即線性系統的維 數), 此時如果直接採用上一節的直接法是一 很不明智的做法, 因為如果這個線性系統的 維數是 n 且沒有特殊的結構, 則需有足夠的 記憶體來儲存約 n2 個實數及約需 n3 的運 算才能得到線性系統的解, 舉例來說, 如果 n = 10000 則約需 4 天的運算時間及 800 mega byte(以一個實數需 8 bytes 來計算)
的記憶體。 到底要如何解這類的問題呢? 我 們之所以用“這類”的字眼, 在於它們通常都 有一共同的特性—–稀疏 (sparsity), 亦即這 類的線性系統的係數所構成的矩陣 A 中大部 分元素都是零, 因此我們可以避免乘零或加 減零的運算來減少運算時間及只存非零的元 素以減少所需的記憶體量。 我們可以直接法 配合上述的技巧來解此類的線性系統, 但是 這需要非常小心及有技巧的處理矩陣 A 的元 素, 因為高斯消去法會把 A 中原來為零的元 素變成非零元素 (亦即 fill-in 的問題), 所以 必須非常技巧的利用 pivoting 來減少 fill-in 的量; 另外由於 fill-in 的問題, 而很難預估記 憶體的需要量, 因此直接法通常並不被建議 來處理大型且稀疏的線性系統, 也因此而發 展出可輕易地避免無謂的乘零及加零的運算 及使用最少記憶體量的迭代法。
根據一般說法, 約 170 年前 Carl Friedrich Gauss 是最早使用疊代法來解 線性系統的人, 那時候他在處理最小方差問 題時, 發現所導出的線性系統太大而無法以 直接法中的高斯消法去求解, 他那時候所用 的方法就是現今大家所熟悉的 Blockwise Gauss-Seidel 疊代。 此後不久, 約1840 Carl Gustav Jacobi 也提出一非常類似的方法, 後來被稱為 Jacobi 疊代。 由於電腦的快速發 展, 人類可以處理更大的線性系統後, 這時才 有人發現及證明上述的兩個方法的收斂速度 太慢, 所以才開始思考如何增加它的收斂速 度, 其中最有名的就是 1950, D. Young 所 提出的 SOR(successive over relaxation) 疊代, 從那時候起, 陸陸續續有人從事這方面 的研究, 提出更多有效且快速的方法。
對一線性系統 Ax = b, 首先選取任意 初始向量 (initial vector) x(0), 利用 x(0) 推 得 x(1), 再以 x(1) 推出 x(2), 依此類推而得 到一向量數列 {x(m)}∞m=0, 這就是迭代法的 主要觀念。 隨著從 x(m) 推演得 x(m+1) 的方 式的不同, 而有形形色色不同的迭代法, 如前 述的 Jacobi 疊代的方式為
x(m+1)i
= (bi−
i−1
X
j=1
aijx(m)j −
n
X
j=i+1
aijx(m)j )/aii
i = 1, 2, . . . , n (7)
Gauss-Seidel 疊代的方式 x(m+1)i
= (bi−
i−1
X
j=1
aijx(m+1)j −
n
X
j=i+1
aijx(m)j
/aii
i = 1, 2, . . . , n (8)
這兩個迭代方法的架構大同小異, 它們都源 於下列的觀念: 每一矩陣 A 都可以唯一分 解成 A = D − L − U 其中 D, L, U 分別為對角, 嚴格下三角及嚴格上三角矩陣, 且 Ax = b, Dx = b + (L + U)x 及 (D − L)x = b + Ux 三個線性系統是相等 的 (equivalent), 所以把左邊及右邊向量 x 分別以 x(m+1) 及 x(m) 代替之, 則可分別得 到上述的 Jacobi 迭代
Dx(m+1) = b + (L + U)x(m) (9) 及 Gauss-Seidel 迭代
(D − L)x(m+1) = b + Ux(m) (10)
從上述的式子, 我們可輕易地發現, 如果 矩陣 D 是奇異的 (singular) 則這兩個迭代 法都無法執行, 然而以上的觀念可推廣為改 考慮矩陣 A 的任意分解 (splitting) A = M − N 而得到下列疊代式
Mx(m+1) = Nx(m)+ b (11) 為了能有效及輕易地執行上述迭代式, 通常 我們會要求 M 是可逆的, 且它的反矩陣 M−1 不難求得(例如: 如果 D 是可逆的, 則 可取 M = D 或 M = D − L 亦即得到所 謂的 Jacobi 及 Gauss-seidel 迭代式)。 在 M 是可逆的情況下, 假設 x∗ 是 Ax = b 的 正確解則由 Mx∗ = Nx∗+ b 我們得到
x(m+1)− x∗
= M−1N(x(m)− x∗) (12)
= (M−1N)2(x(m−1)− x∗)
= · · ·
= (M−1N)m+1(x(0)− x∗)
從式子 (12) 及利用線性代數的定理, 我們 可歸納出, 如果 M−1N 的 spectral 半 徑 ρ(M−1N) < 1, 則 {x(m)}∞m=0 是 一收斂數列, 且其收斂向量 ˜x = x∗, 如 果 ρ(M−1N) > 1 則 {x(m)}∞m=0 必 為一發散數列, 如果 ρ(M−1N) = 1 則 {x(m)}∞m=0 可能收斂也可能發散。 另外我們 亦可以 ρ(M−1N) 的大小來粗略的估計其 收斂速度, 前述所提 Jacobi 或 Gauss- Seidel 迭代的收斂太慢的原因就是它們相對 應的迭代矩陣 M−1N 的 spectral 半徑
ρ(D−1(L + U)) 及 ρ((D − L)−1U) 不夠 小。 為了使其收斂快些, 我們引進參數 ω, 而 改考慮線性系統 ωAx = ωb, 且考慮分解 M = D − ωL 及 N = (1 − ω)D + ωU 而 得到所謂的 SOR 迭代式
(D−ωL)x(m+1) = ωb+[(1−ω)D+ωU]x(m) (13) 亦即
x(m+1)i
= x(m)i −ω
i−1
X
j=1
aijx(m+1)j +
n
X
j=i
aijx(m)j −bi
/aiii = 1, 2, · · · , n (14)
因為 ω 是一任意實數, 且 ω = 1 即所謂的 Gauss-Seidel 迭代, 所以一定存在 ω 使得 ρ((D−ωL)−1((1−ω)D+ωU)) ≤ ρ((D−
L)−1U), 其中使得 ρ((D − ωL)−1((1 − ω)D + ωU)) 為最小的 ω, 我們稱之為最佳 的 ω, 以 ωopt 代表之。 找 ωopt 並非一簡單 的工作, 在1950及1960年代有不少的學者從 事這方面的研究, 於某些特定類型的矩陣找 到了 ωopt, 我們不在此多作描述, 有興趣的 讀者可參照後面所列的參考資料。
四 . 結語
進入數值線性代數的領域不難, 只要具 備線性代數的基本概念及些許的程式設計能 力即可, 一旦你登堂入室後, 你便可發現其內 部的奧妙, 它並不如想像中的容易, 因為還是 有很多看似淺顯的問題尚未被解決出, 留待 聰明的讀者們去解決。
在此我們僅介紹兩個淺顯易懂的方法給 讀者們, 希望能引起讀者的興趣而動手寫些 程式去執行上述的兩大演算法。 唯有動手做 問題, 您才會發現其實還有很多的問題我們 在此尚未探討的, 諸如直接法中的捨入誤差 對解的影響有多大, 迭代法中須如何判斷數 列 {x(m)}∞m=0 是否收斂到正確解了, 及怎 樣的判斷得到的解最接近正確解等等的問題。
希望做更深入的探討的讀者, 不妨參考我們 所列的參考資料。
