在科學計算和工程應用領域，經常涉及求解線性方程組問題。迭代法是解決這類問題的有力武器，通過多次迭代逐步逼近精確解，在處理大型稀疏線性方程組時具有重要優勢。

撰文|丁玖（美國南密西西比大學數學系教授）

記得有個大科學家曾經說過，這個世界是非線性的。然而非線性方程一般不能直接求解，即解析解雖可證存在卻無具體表達式，因而迭代法幾乎是唯一可行的辦法。比如說，從介值定理可知，方程x = cos x在區間(0, 1)內定有一解，但沒有一步到位的法子找到它，人們只能用基於介值定理的二分法或基於切線逼近的牛頓法，來求得此方程的迭代近似解。這樣，從最古老的巴比倫平方根迭代法，到今日非線性方程組數值解的最重要方法——牛頓迭代法，人們一直熱衷於迭代法的理論探索和創新實踐。

既然有非線性方程，就必然有線性方程，而且後者在科學技術和日常生活中同樣隨處可見。電話公司的客戶安排、航空公司的航班分配、運輸公司的運貨路線……都可以用線性規劃來最佳化成本。四十年前，當我在南京大學數學系時，讀到一篇關於線性規劃的綜述性文章，其中寫道：全世界的計算機大約有百分之六十的時間，都在求解線性規劃問題！

另一方面，微分概念告訴我們，非線性函數可以局部地用線性函數來近似，構造牛頓迭代法本質上就基於這個簡單的觀察。只要非線性函數f在一個點a存在導數，那麼在a點的附近，函數值f(x)就約等於線性函數值f(a) + f'(a)(x-a)，其幾何直觀則是函數圖像的切線在切點附近與曲線相差無幾。

此外，在物理科學中，許多描述自然規律的常微分方程或偏微分方程對於未知函數以及它們的導數或偏導數都呈現出線性關係，如熱傳導方程或波動方程。當用計算數學中的差商逼近導數技術或有限元思想來離散化這些連續方程時，獲得的差分方程也是線性的。它們可以寫成一般形式Ax = b，其中A是一個n行、n列的係數矩陣，b是一個n維的已知列向量，x是一個n維的未知列向量。這個線性方程看上去像一元一次方程ax = b一樣簡單，但如果按照矩陣乘法的法則將方程左邊每個分量的代數表達式全部寫出來，結果就是一組含有n個未知數x₁, x₂, …, x_n的n個n元一次方程。如果將方矩陣A中第i行、第j列的元素記為a_ij，將列向量b的第i個分量記為b_i，那麼線性方程組Ax = b展開後的第i個方程為

a_i1x₁+ a_i2x₂+ … + a_inx_n= b_i, i = 1, 2, …, n。

計算數學中的一個子學科「數值線性代數」有個基本的職責：怎樣有效求解上面的線性方程組？當n = 2或3，最多不超過4時，我們早已在中學課堂上學會了怎樣求解二元一次方程組、三元一次方程組，等等。採用的方法無非是兩種，一種叫做代入法，另一種稱為加減法。它們都是基於同一個思想：通過消去方程中的一個或幾個未知數，設法獲得只含有一個未知數的方程，比方說-2y = 3，從而解出這個未知數。

中學生所學的上述方法，實際上是大名鼎鼎的高斯消去法的簡單實現，發明者正是德國天才高斯 (Carl Friedrich Gauss，1777-1855) 。高斯消去法的基本思想還是如上所說：將含有多個未知數的方程化簡到只剩下一個未知數，因為只含一個變元的線性方程如-2y = 3總是可以一步到位解出y = -3/2的。基於這個非常簡單的思路，高斯想出了一個點子，將線性方程組的係數方陣程式化系統性地化約成主對角線下方的元素全為零的上三角矩陣，這樣新的等價的線性方程組的最後一個方程只含一個未知數，倒數第二個方程含有兩個未知數，如此等等。因而，通過所謂的「回代法」，從最後一個方程解出唯一的未知數，然後將它的值代入到倒數第二個方程，結果該方程也只含一個未知數，馬上可以解出，再將這兩個未知數的值代入到倒數第三個方程，同法解出下一個的未知數。如此這般，就可以依次解出所有的未知數。

高斯消去法在經過有限次代數運算後就可以完全找出給定線性方程組的精確解（如果每步運算的誤差為零的話）。所以它屬於求解線性方程組的第一類數值方法——直接法，這是解線性方程組優越於解非線性方程組的一個地方。然而，對有巨量（如成千上萬個或更大數量級）變元個數或者特別多的係數為零（比如係數矩陣為三對角或其他稀疏類型）的那種方程組，直接法常常沒有第二類方法——迭代法有效。對線性偏微分方程直接用差商代替偏導數的差分方法，或將微分方程問題變為數學上等價的變分問題而進行數值最佳化的有限元方法，所得到的大型線性方程組一般都帶有特定結構的稀疏係數矩陣，這時，迭代法幾乎是最有效的演算法了。

我們下面只講解線性方程組的迭代法。讓我們回憶求解單變數非線性方程的迭代法，一般形式是x_n= f(x_n-1), n = 1, 2, …，其中f是將定義域區間映到自身的一個函數，x₀是迭代所取的初始點。然而，對於線性迭代法，迭代函數不再是一個自變數的線性函數，而是有n個自變數的線性向量函數。由於字母n現在另有他用，我們將用字母k代表迭代次數的下標，而將多變數線性向量函數用y = L(x)表示，其中L(x)的表達式是Mx + c，M是一個有n行和n列的給定矩陣（也稱為n階正方矩陣或n階方陣），c是一個給定的n維列向量，其分量是c₁, c₂, …, c_n，x = (x₁, x₂, …, x_n)^T是n維變元列向量，其中的上標字母T表示轉置運算，即一般的有m行和n列的矩陣A的轉置矩陣A^T有n行和m列，其第j行的元素是原矩陣第j列的對應元素。

本章將要面對的線性迭代程序可以寫成

x_k= L(x_k-1) = Mx_k-1+ c, k = 1, 2, …; x₀是給定的。

我在以前的科普文章裡寫過，單變數非線性方程的迭代法收斂的一個充分條件，是迭代函數在不動點的導數絕對值小於1。那麼，如果被迭代的是有多個變數的線性向量函數，什麼樣的條件可以充分到保證迭代法收斂呢？所謂序列的收斂，本質上是用距離這個概念來定義的，這個概念我們在以前提到巴拿赫不動點定理或壓縮對映定理時已經見到過。我們再來回憶有關內容。在一個抽象的距離空間(X, d)裡，一個由X中的元素所組成的無窮序列{x_k}被說成是收斂到X中的一個元素x，指的是由距離函數d所確定的非負數列{d(x_k, x)}當k趨向於無窮大時趨向於0。對於單變數函數迭代法的收斂性問題，我們實際上選取了一個簡單而具體的距離空間，它就是所有實數全體構成的一維歐幾里得空間R，其任意兩個實數之間的距離就是實數軸上兩點x和y之間的通常距離|x – y|。

上面兩數之間的距離可用線段長度的說法來等價地定義。任意一個實數a的絕對值|a|在幾何上的意義是以數軸上代表a的那個點和原點0為兩端點的線段的長度。不難看出，兩個數x和y之間的距離就是數x – y所對應的那條線段的長度。

對於現在要考察的高維歐幾里得空間Rⁿ，兩個n維向量之間的距離可如法炮製地定義。對於二維平面R²，設第一個行向量為x = (a, b)，第二個行向量是y = (c, d)，則它們之間的距離就是平面上這兩點之間的通常距離

， 或等價地，它是平面向量x – y的長度。在通常的三維空間R³裡，設向量x 的三個分量是a, b, c，向量y的三個分量是u, v, w，則x和y之間的距離也為這兩個空間點之間的通常距離

。

它也就是空間向量x – y的實際長度。

三維以上的歐幾里得空間，我們的眼力無法感知，但我們的想象力可以穿越界限，所以將上面直觀的距離公式推廣開來，就得到n維空間Rⁿ的歐幾里得距離定義：設向量x的分量為a₁, a₂, …, a_n，向量y的分量為b₁, b₂, …, b_n，則它們之間的距離是

d(x, y) = {(a₁– b₁)²+ … + (a_n– b_n)²}^1/2。

我們將三維以上空間中的向量「長度」改稱為「歐幾里得範數」，簡稱範數或2-範數，則n維向量x的範數||x||₂被定義為x的所有分量的平方之和再開平方根。這樣，x和y之間的距離就等於向量x – y的歐幾里得範數||x – y||₂。

作為長度概念的推廣，範數繼承了長度的三個基本性質：

1. 向量x的範數||x||₂總是非負實數，且||x||₂= 0當且僅當x是零向量，即它的每個分量都為0。

2. 一個數α乘以向量x得到的向量αx，其範數等於α的絕對值乘以x的範數，即||αx||₂= |α| ||x||₂。這裡αx的各個分量等於α乘以x的相應分量。

3. 向量x和向量y的和的範數小於或等於向量x的範數加上向量y的範數，即||x + y||₂≤ ||x||₂+ ||y||₂。這裡x + y的每個分量為x和y的對應分量之和。

性質(iii)中的不等式一般被稱為三角不等式，因為在二維和三維的情形，它就是平面幾何中的一條定理：三角形的任意一邊的長度不大於其餘兩邊的長度之和。

類似於關於絕對值的三角不等式，範數的三角不等式可以匯出不等式| ||x||₂– ||y||₂| ≤ ||x – y||₂。由此可知將x映到 ||x||₂的範數函數|| ||₂是定義域Rⁿ上的一個連續函數。

有了n維歐幾里得空間Rⁿ的向量範數概念，定義在Rⁿ上的多變數線性向量函數L可以按照線性代數或泛函分析的習慣叫法被稱為仿射運算元；特別地，當常向量c = 0時，L叫做線性運算元。這裡，「線性」意味著對所有的數α和β以及向量x和y，恆有等式L(αx + βy) = αL(x) + βL(y)，而「仿射」則意味著對所有的滿足α + β = 1的數α和β以及任意向量x和y，恆有等式L(αx + βy) = αL(x) + βL(y)。這樣，我們可以用上由向量範數所誘匯出的矩陣範數，結合著名的「壓縮對映定理」來檢驗線性迭代法的收斂性。為此目的，我們將Rⁿ中所有其歐幾里得範數不大於1的那些向量組成的集合稱為Rⁿ的n維單位閉球，它的「表面」，即所有範數等於1的向量全體稱為(n – 1)維單位球面，用S^n-1表示。

給定一個n階方陣M，讓x取遍Rⁿ的(n – 1)維單位球面中的所有向量，則定義了一個非負實函數，它把x映到||Mx||₂。由於線性運算x → Mx和範數運算u → ||u||₂都是連續的，它們的複合函數||Mx||₂在定義域S^n-1上是處處連續的。又因為有限維的單位球面不僅是有界集合（因為集合中所有元素的範數一致有界；1就是一個上界），而且是閉集（因為只要S^n-1中的任意一個向量序列{x_k}當k趨向於無窮大時收斂到向量x，則在恆等式||x_k||₂≡ 1中兩邊取k趨向於無窮大時的極限，考慮到範數的連續性，就有||x||₂= 1，即x也屬於S^n-1），因此根據微積分中的一條定理：定義在歐幾里得空間的有界閉集上的連續函數一定有最大值，我們得出結論：非負函數||Mx||₂在S^n-1上的某一點取到最大值。沒有運氣學過連續函數最大值定理的讀者肯定在中學學過代數，可以通過畫出站在有界閉區間[-1, 1]上的一條包含兩個端點的拋物線，馬上感知：此連續曲線上有一點最高（當然也有一點最低）。

這個最大值就被定義為由Rⁿ上的歐幾里得範數所誘匯出的矩陣範數||M||₂。從矩陣運算的性質馬上得知如下關於矩陣、向量範數的常用不等式

||Mx||₂≤ ||M||₂||x||₂對所有的n維列向量x都成立。

另外，被向量範數誘匯出的矩陣範數也和向量範數一樣滿足上述三項「範數」的基本性質，即||M||₂是非負實數，它的值為零當且僅當M是零矩陣；一個數α和M的數乘αM的範數等於α的絕對值乘以M的範數；兩個同階矩陣之和的範數不大於它們各自的範數之和。由於兩個同階方陣可以相乘，矩陣範數比向量範數多了一個基本性質：兩個同階方陣之積的範數不大於它們各自的範數之積。特別地，||M^k||₂≤ (||M||₂)^k。

是時候尋找巴拿赫不動點定理的新應用了。對於線性迭代向量函數L(x) = Mx + c，我們有

||L(x) – L(y)||₂= ||Mx + c – (My + c)||₂= ||M(x – y)||₂≤ ||M||₂||x – y||₂。

如果矩陣M滿足嚴格不等式

||M||₂ 1，那麼，L在完備的距離空間(Rⁿ, d)上是壓縮的；這裡的距離d(x, y) = ||x – y||₂。空間Rⁿ的「完備性」是由實數集R的完備性繼承而來。根據巴拿赫不動點定理，方程x = L(x)有唯一的不動點列向量p，並且從任一個初始列向量x₀出發，迭代點列x_k= L(x_k-1), k = 1, 2, …當k趨向於無窮大時將收斂到p。這就是線性迭代法的一個收斂性條件。自然，這僅僅是一個充分條件，不是必要條件。在後續的一篇文章裡，我們將想方設法推匯出線性迭代收斂的一個既充分又必要的條件。

我們將如上的充分條件用到迭代求解線性方程組Ax = b，其中的方陣A假定為非奇異的，這保證了對任意右端列向量b，該方程組有唯一解。先解釋什麼叫非奇異矩陣。一個有n行和n列的矩陣M可以被看成是定義在n維歐幾里得空間Rⁿ上的一個線性運算元，它將每一個n維列向量x映到n維列向量Mx。自然它將n維零向量0映到它自己（故向量0是M的一個不動點）。如果M不能將一個非零列向量（即該向量至少有一個分量為非零數）映成向量0，或等價地說，如果M將所有的非零向量映成非零的向量，則我們稱M為一個非奇異矩陣。非奇異矩陣又被稱為可逆矩陣，因為一個方陣是非奇異的當且僅當它一一對應地將Rⁿ中的每一個元素映到Rⁿ中的每一個元素。換言之，可逆矩陣不僅是單射的（即不同的向量被映到不同的向量），而且是滿射的（即Rⁿ中的每一個列向量都是該矩陣乘以某一個列向量的結果）。

如果n階方陣M是非奇異的，那麼它是單射和滿射的。單射和滿射這兩個性質同時成立時就叫做是雙射的，因而M作為定義域為Rⁿ的線性運算元，它具有所謂的逆運算元，其定義域是M的值域Rⁿ。故存在一個與M同階的方陣，稱為M的逆矩陣，記為M^-1，滿足條件：對任意的n維列向量x和y，等式Mx = y成立當且僅當等式M^-1y = x成立。如果把運算元M視為一個函數，那麼逆矩陣實質上無異於中學代數課本里講到的反函數概念。如此一來，可逆矩陣和它的逆矩陣同時滿足兩個等式MM^-1= I及M^-1M = I，這和可逆函數f及其反函數f^-1滿足的兩個恆等式f^-1(f(x)) ≡ x，f(f^-1(y)) ≡ y一模一樣，其中x取f的定義域中的所有元素，y取f的值域中的所有元素。這些都是最簡單的算術等式aa^-1= a^-1a = 1(a ≠ 0)的自然推廣。我們也順便指出，從矩陣理論可知，上述關於矩陣和逆矩陣的兩個等式並非彼此獨立，其中的一個成立就能推出另一個也成立。如果讀者熟悉行列式的概念，可能也會知道方陣為非奇異的一個簡單的等價條件：矩陣的行列式不等於0。

有了上面的準備知識，我們可以回到設計迭代法求解具有標準形式的線性代數方程組Ax = b這個文章主題。第一個需要考慮的問題是怎樣將此向量方程寫成迭代法的標準形式x = L(x)。方法很簡單，只要將矩陣A分解為兩個矩陣之差就行：A = N – P，但需要一個額外條件，那就是矩陣N必須是非奇異的。然後原方程Ax = b就可寫成Nx = Px + b。由於N是可逆矩陣，前述方程等價於下列不動點方程

x = N^-1Px + N^-1b = Mx + c，其中M = N^-1P和c = N^-1b。因此，對應的求解線性方程組Ax = b的迭代公式是

x_k= N^-1Px_k-1+ N^-1b, k = 1, 2, 3, …; 給定x₀。

讀者自然會問，怎麼選取關鍵的矩陣N？這個矩陣要有逆矩陣，要保證迭代法收斂，並且是收斂得越快越好。這是計算數學家們一直追求的目標。我們現在潛入歷史，看看有哪幾位數學家線上性迭代法方面貢獻卓越，成為先驅。

歷史上第一個提出迭代法的數學家又是高斯，所以他不僅是以高斯消去法為代表的直接法的鼻祖，也是迭代法的鼻祖。1823年12月26日，高斯在給他昔日博士生戈林 (Christian Ludwig Gerling，1788-1864) 的一封信中，提出了這個他稱之為「非直接消去法」的迭代法，用於求解關於四個城市的一個應用問題。他首先利用也是他發明的最小二乘法得到一組線性方程，然後用自己設計出的迭代法求解對應的法方程。這走出了迭代法求解線性方程組歷史進程的第一步。然而，由於高斯的許多數學發現和發明並沒有被他正式發表，所以他的發明權不得不常與後人分享，他也不在乎，沒有當今「不發表便滅亡」的「科學家煩惱」。半個世紀後，他的同胞數學家馮·賽德爾 (Phillip Ludwig von Seidel，1821-1896) 於1874年發表了包含同樣思想的迭代法。

1845年，另一名德國數學家雅可比 (Carl Gustov Jacob Jacobi，1804-1851) 也提出了一個以他的名字命名的迭代法。這個方法比較簡單，在數值代數的教科書中一般是作為迭代法求解線性方程組的第一個例子被介紹給學生。雅可比的人生軌跡雖然只持續了不到四十七年，他卻是一位多產的數學大師，其最偉大的貢獻在橢圓函數方面。他被上世紀的美國數學史家貝爾 (Eric Temple Bell，1883-1960) 寫進影響了幾代讀者的《大數學家》 (Men of Mathematics)，該書總共只描述了到上世紀初為止留名數學史的三十多個「大數學家」。雅可比的職業生涯中有十三年是在柯尼斯堡大學擔任正教授，這所大學的傑出校友有數學家希爾伯特 (David Hilbert，1862-1943) 和閔可夫斯基 (Hermann Minkowski，1864-1909)。

求解線性方程組歷史上最有名氣也最簡單實用的兩個迭代法都是德國人的傑作，今日關於它們的基本理論也很漂亮。但是本文已夠長了，所以我將再寫一文，專門介紹雅可比迭代法和高斯-賽德爾迭代法的基本思想和收斂性。

寫於2023年11月3日星期五

美國哈蒂斯堡夏日山莊

本文受科普中國·星空計劃項目扶持

出品：中國科協科普部

監製：中國科學技術出版社有限公司、北京中科星河文化傳媒有限公司

本文經授權轉載自微信公眾號：返樸 作者：丁玖