4 線性方程組的直接解法
4.1 引言
在自然科學和工程計算的很多問題通常都可以歸結為求解線性方程組的問題，如三次樣條插值、最小二乘法擬合曲線等等。因此，在本章中將探討線性方程組的直接解法。
一般的n元線性方程組具有以下形式：
\[\left\{\begin{aligned} a_{11}x_1+a_{12}x_2&+\cdots+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2&+\cdots+a_{2n}x_n=b_2\\ &\cdots\\ a_{n1}x_1+a_{n2}x_2&+\cdots+a_{nn}x_n=b_n \end{aligned}\right. \]
寫成矩陣的形式就是：\(Ax=b\)
其中
\[A=\left[\begin{matrix} a_{11}&a_{12}&\cdots&a_{1n}\\ a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n1}&a_{n2}&\cdots&a_{nn}\\ \end{matrix}\right], x=\left[\begin{matrix} x_1\\x_2\\\vdots\\x_n \end{matrix}\right] ,b=\left[\begin{matrix} b_1\\b_2\\\vdots\\b_n \end{matrix}\right] \]
如果線性方程組的系數行列式不為零，即\(\det(A)\ne0\)，那么該方程有唯一解。根據克萊姆（Gramer）法則，方程的解為
\[x_i=\frac{\det{A_i}}{\det{A}},i=1,2,\cdots,n \]
如果使用克萊姆法則求解，那么一共需要計算\(n+1\)個\(n\)階行列式，並進行\(n\)次除法。計算一個\(n\)階行列式需要\((n-1)\times n!\)次乘法，總的計算量就是\(n+(n-1)(n+1)!\)，當\(n\)的值比較大的時候，需要執行的浮點數計算次數過於巨大，因此在實際中是不可行的。
關於線性方程組的數值解法一般有兩類：
- 直接法：經過有限步運算，求得方程組的精確解（如果不考慮計算過程中的舍入誤差）
- 迭代法：用某種極限過程去逼近線性方程組的精確解，程序設計比較簡單，但存在收斂性和收斂速度等問題。
4.2 消元法
4.2.1 三角形方程組
首先，介紹兩種最簡單方程組形式：上三角形方程組、下三角形方程組。實際上消元法的目的就是將一般的方程組轉化為三角形方程組進行求解。
上三角形方程組的一般形式為：
\[\left\{\begin{aligned} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n&=b_1\\ a_{22}x_2+\cdots+a_{2n}x_n&=b_2\\ \vdots&\\ a_{nn}x_n&=b_n \end{aligned}\right. \]
其中\(a_{ii}\ne0(i=1,2,\cdots,n)\)。
為了求解上三角形方程組，我們先求解出\(x_n=b_n/a_{nn}\)，然后就可以從下至上利用已經求出的\(x_k\)，依次求出\(x_{n-1},\cdots,x_1\)的值，這就完成了上三角形方程組的求解。整個求解過程可以被描述為：
\[\left\{\begin{aligned} x_n&=\frac{b_n}{a_{nn}}\\ x_i&=\frac{1}{a_{ii}}(b_i-\sum_{k=i+1}^n{a_{ik}x_k}), (i=n-1,\cdots,2,1) \end{aligned}\right. \]
上述求解過程稱為回代，因此也將這種方法稱為回代法。整個過程需要進行乘法次數計算如下
\[1+2+\cdots+(n-2)+(n-1)=\frac{n(n-1)}2 \]
需要執行的除法次數為\(n\)次，那么一共就需要\(n(n+1)/2\)次浮點運算。
同理，對於下三角形方程組的求解方法為
\[\left\{\begin{aligned} x_1&=\frac{b_1}{a_{11}}\\ x_i&=\frac{1}{a_{ii}}(b_i-\sum_{k=1}^{i-1}{a_{ik}x_k}), (i=2,3,\cdots,n-1) \end{aligned}\right. \]
所需的浮點運算次數和求解上三角形方程組相同。
4.2.2 高斯消元法
消元法的基本思想就是對方程組進行初等變換，把一般形式的方程組轉化為容易求解的三角形方程組。基本過程就是先將一般的線性方程組轉換成通解的上三角形方程組（這個過程稱為消元），然后利用回代法即可求解原方程組的解。
消元過程如下所示：
假設現有一個\(n\)階方程組（\(a_{11}\ne0\)）：
\[\left\{\begin{aligned} a_{11}x_1+a_{12}x_2&+\cdots+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2&+\cdots+a_{2n}x_n=b_2\\ &\cdots\\ a_{n1}x_1+a_{n2}x_2&+\cdots+a_{nn}x_n=b_n \end{aligned}\right. \]
令\(m_{i1}=-a_{i1}/a_{11}\)，用\(m_{i1}\)乘以第一個方程，然后加到第\(i\)個方程上，得到
\[\left\{\begin{aligned} a_{11}x_1+a_{12}x_2&+\cdots+a_{1n}x_n=b_1\\ a_{22}^{(1)}x_2&+\cdots+a_{2n}^{(1)}x_n=b_2^{(1)}\\ &\cdots\\ a_{n2}^{(1)}x_2&+\cdots+a_{nn}^{(1)}x_n=b_n^{(1)} \end{aligned}\right. \]
其中\(a_{ij}^{(1)}=a_{ij}+a_{1j}m_{i1},b_i^{(1)}=b_i+b_1m_{i1}(i,j=2,3,\cdots,n)\)。
如果\(a_{22}^{(1)}\ne0\)，則保留前兩個方程。與上一步類似，進一步進行消元，得到
\[\left\{\begin{aligned} a_{11}x_1+a_{12}x_2+a_{13}x_3+\cdots+a_{1n}x_n&=b_1\\ a_{22}^{(1)}x_2+a_{23}^{(1)}x_3+\cdots+a_{2n}^{(1)}x_n&=b_2^{(1)}\\ a_{33}^{(2)}x_3+\cdots+a_{3n}^{(2)}x_n&=b_3^{(2)}\\ \vdots&\\ a_{n3}^{(2)}x_2+\cdots+a_{nn}^{(2)}x_n&=b_n^{(2)} \end{aligned}\right. \]
其中\(a_{ij}^{(1)}=a_{ij}^{(1)}+a_{2j}m_{i2},b_i^{(2)}=b_i^{(1)}+b_2^{(1)}m_{i2}(i,j=3,4,\cdots,n),m_{i2}=-a^{(1)}_{i2}/a_{22}^{(1)}\)。
不斷重復上述過程，最終可以得到等價的上三角方程組：
\[\left\{\begin{aligned} a_{11}x_1+a_{12}x_2+a_{13}x_3+\cdots+a_{1n}x_n&=b_1\\ a_{22}^{(1)}x_2+a_{23}^{(1)}x_3+\cdots+a_{2n}^{(1)}x_n&=b_2^{(1)}\\ a_{33}^{(2)}x_3+\cdots+a_{3n}^{(2)}x_n&=b_3^{(2)}\\ \vdots&\\ a_{nn}^{(n-1)}x_n&=b_n^{(n-1)} \end{aligned}\right. \]
此時就可以使用回代法進行求解。
現在我們來考慮進行消元所需要執行的乘除操作次數。考慮第\(k\)步消元時所需要的操作，首先，計算\(m_{ik}(i=k+1,\cdots,n)\)需要\(n-k\)次除法。然后，將第\(k\)個方程乘以\(m_{ik}\)需要\(n-k+1\)次乘法（不是\(n-k+2\)，因為\(a_{kk}^{(k-1)}\)下方的數均是0，可以直接賦值不需要實際計算）。因為一共有\(n-k\)個方程，所以最終所需的浮點操作為\((n-k)+(n-k)(n-k+1)\)。消元步驟一共執行\(n-1\)次，所以全部消元過程所需的乘除操作次數為
\[\sum_{k=1}^{n-1}=[(n-k)+(n-k)(n-k+1)]=\frac{1}{6}n(n-1)(2n+5) \]
而回代過程所需的乘除操作次數為\(n(n+1)/2\)，所以高斯消元法所需的總計算量為
\[N=\frac{1}{6}n(n-1)(2n+5)+\frac{n(n+1)}{2}=\frac{1}{3}n(n^2+3n-1) \]
因此計算復雜度為\(O(n^3)\)，所以采用高斯消元法的計算量將遠小於使用克萊姆法則所需的計算量。
這里還需要主要的問題是，在上述高斯消元的過程中，我們是從上至下順序進行的，這就要求在\(a_{kk}^{(k-1)}\)處的值不等於0，即對系數矩陣主對角線上的元素不為零。換句話說，如果高斯消元法能夠順利進行就需要保證：\(a_{kk}^{(k-1)}\ne0,k=1,2,\cdots,n\)。
考慮系數矩陣\(A\)的順序主子矩陣的行列式：
\[\det A_k=\left|\begin{matrix} a_{11}&a_{12}&\cdots&a_{1k}\\ a_{21}&a_{22}&\cdots&a_{2k}\\ \vdots&\vdots&\ddots&\vdots\\ a_{k1}&a_{k2}&\cdots&a_{kk}\\ \end{matrix}\right|= \left|\begin{matrix} a_{11}&a_{12}&\cdots&a_{1k}\\ 0&a_{22}^{(1)}&\cdots&a_{2k}^{(1)}\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&a_{kk}^{(k-1)}\\ \end{matrix}\right|= a_{11}a_{22}^{(1)}\cdots a_{kk}^{(k-1)} \]
其中\(k=1,2,\cdots,n\)。
如果\(A\)所有的順序主子式\(\det A_k\ne0\)，就可以保證在\(a_{kk}^{(k-1)}\)處的值不等於0，反之亦然。因此順序高斯消元法可行的充要條件為\(A\)的所有順序主子式均不為0。
4.2.3 列主元消元
上述高斯消元法存在一定的問題，下面給出兩個例子：
例題1 使用高斯消元法求解如下方程組：
\[\left\{\begin{aligned} 2x_2=1\\ 2x_1+3x_2=2 \end{aligned}\right. \]
此時\(a_{11}=0\)，不能直接使用高斯消元法。需要先交換兩個方程的位置才行。
例題2 假設我們在，只能表示四位浮點十進制數計算機上計算下列方程組
\[\left\{\begin{aligned} 0.0001x_1+x_2=1\\ x_1+x_2=2 \end{aligned}\right. \]
本題在計算機內部表示為
\[\left\{\begin{aligned} 0.1000\times10^{-3}x_1+0.1000\times10^1x_2&=0.1000\times10^1\\ 0.1000\times10^1x_1+0.1000\times10^1x_2&=0.2000\times10^1 \end{aligned}\right. \]
此時是可以直接使用高斯消元法的，因為\(a_{11}\ne0,m_{21}=-10^4\)，所以
\[\begin{aligned} a_{22}^{(1)}&=0.1000\times10^1-10^4\times0.1000\times10^1\\ &=0.00001\times10^5-0.1000\times10^5\\ &=0.0000-0.1000\times10^5\\ &=-0.1000\times10^5 \end{aligned} \]
由於浮點數計算時，需要進行對階，這使得小數被大數“吃掉了”。同理我們可以得到
\[b_2^{(1)}=-0.1000\times10^5 \]
於是得到三角方程組
\[\left\{\begin{aligned} 0.1000\times10^{-3}x_1+0.1000\times10^1x_2&=0.1000\times10^1\\ -0.1000\times10^5x_2&=-0.1000\times10^5 \end{aligned}\right. \]
通過回代可以解得\(x_2=1,x_1=0\)。
然而原方程的實際解為\(x_2=9998/9999,x_1=10000/9999\)，因此此題如果直接使用高斯消元法將會導致解嚴重失真。導致這個問題出現的原因在於消元過程中除了一個很小的數使得最后的結果非常大，這可能導致出現大數吃小數的現象，如果原問題是病態的將會導致最終結果嚴重偏離理想值。
為了解決這個問題，將引入選主元消元法，避免除以一個很小的數。在高斯消元過程中，我們將\(a_{kk}^{(k-1)}\)位置稱為主元位置，位於該位置的元素則稱為主元素。如果在進行第\(k\)步消元時，將第\(k\)列中絕對值最大的元素作為主元素，將其調換到第\(k\)行在做消元，這種方法就稱為列主元消元。具體方法如下：
- 先在第一列選取絕對值最大的元素：\(|a_{m1}|=\max\{|a_{11}|,|a_{21}|,\cdots,|a_{n1}|\}\)。交換第1行和第m行的所有元素（即交換第1個方程和第m個方程的位置），在做消元。
- 進行第\(k(\ge2)\)步消元前，在選取第\(k\)列絕對值最大的元素：\(|a_{mk}^{(k-1)}|=\max\{|a_{kk}^{(k-1)}|,|a_{(k+1)k}^{(k-1)}|,\cdots,|a_{nk}^{(k-1)}|\}\)。交換第1行和第m行的所有元素（即交換第k個方程和第m個方程的位置），在做消元。
與高斯消元法相比，列主元消元只是在每一步消元之前都增加了選主元和交換方程位置的操作，並不會影響最終的結果。
4.2.4 全主元消元
以第\(k\)步消元為例，列主元消元關注的僅僅只是\(a_{kk}^{(k-1)}\)及其下方的元素，而全主元關注的是以\(a_{kk}^{(k-1)}\)作為第一行第一列元素的右下子矩陣的所有元素，因此尋找主元的位置更廣。
如上圖所示，綠色框內表示列主元消元關注的范圍，紅色框內表示全主元關注的范圍。
在執行第\(k\)步消元前，首先，需要在右下子矩陣中選擇絕對值最大的元素作為主元素：
\[|a_{pq}^{(k-1)}|=\max_{k\le i,j\le n}|a_{ij}^{(k-1)}| \]
選定好主元素后，就涉及到上下方向方程位置的互換，以及左右方向未知數位置的互換：
- 若\(p\ne k\)，則交換第\(k\)行和第\(p\)行。
- 若\(q\ne k\)，則交換第\(q\)列和第\(k\)列。
最后，執行消元操作。
注意： 由於列交換改變了未知數\(x_i\)的順序，因此每交換一次，就需要做一次記錄，等求出最終結果之后再交換回來。
4.3 矩陣的分解
4.3.1 LU分解
LU分解的目的是將一個矩陣\(A\)分解成一個下三角陣\(L\)，和上三角陣\(U\)。一旦求得\(A=LU\)，則方程\(Ax=b\)就可以轉化為求解\(LUx=b\)，令\(Ux=y\)，則\(Ly=b\)。實際求解過程如下：
- 由\(Ly=b\)求解出\(y\)。
- 再由\(Ux=y\)求解出\(x\)。
矩陣\(A\)能夠進行LU分解的充要條件是\(A\)的各階主子式不為零，這可以理解為LU分解是通過高斯消元法誘導出來的。
矩陣的LU分解一共有兩種形式：
- 多利特爾（Doolittle）分解：\(L\)為單位下三角陣（主對角線值全為1），\(U\)為上三角陣。
- 庫朗（Crout）分解：\(L\)為下三角陣，\(U\)為單位上三角陣（主對角線值全為1）。
在這里我們只考慮多利特爾分解，下面我們將通過高斯消元法的觀點推導出矩陣的LU分解。
假設\(A\in\mathbb{R}^{3\times3}\)。高斯消元法的每一步消元操作都等價於左乘一個矩陣做初等行變換，首先，考慮第一次消元：
\[L_1=\left[\begin{matrix} 1&0&0\\m_{21}&1&0\\m_{31}&0&1 \end{matrix}\right],m_{i1}=-a_{i1}/a_{11},i=2,3 \]
於是我們有
\[L_1A=\left[\begin{matrix} 1&0&0\\m_{21}&1&0\\m_{31}&0&1 \end{matrix}\right] \left[\begin{matrix} a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33} \end{matrix}\right]= \left[\begin{matrix} a_{11}&a_{12}&a_{13}\\&a_{22}^{(1)}&a_{23}^{(1)}\\&a_{32}^{(1)}&a_{33}^{(1)} \end{matrix}\right] \]
進行第二次消元
\[L_2=\left[\begin{matrix} 1&0&0\\0&1&0\\0&m_{32}&1 \end{matrix}\right],m_{i2}=-a_{i2}/a_{22},i=3 \]
可以得到
\[L_2L_1A= \left[\begin{matrix} 1&0&0\\0&1&0\\0&m_{32}&1 \end{matrix}\right] \left[\begin{matrix} a_{11}&a_{12}&a_{13}\\&a_{22}^{(1)}&a_{23}^{(1)}\\&a_{32}^{(1)}&a_{33}^{(1)} \end{matrix}\right]= \left[\begin{matrix} a_{11}&a_{12}&a_{13}\\&a_{22}^{(1)}&a_{23}^{(1)}\\&&a_{33}^{(2)} \end{matrix}\right]=U \]
根據線性代數中的知識，我們知道初等行變換都是可逆的。因此\(L_1,L_2\)都是可逆矩陣，所以
\[A=L_1^{-1}L_2^{-1}U= \left[\begin{matrix} 1&0&0\\-m_{21}&1&0\\-m_{31}&0&1 \end{matrix}\right] \left[\begin{matrix} 1&0&0\\0&1&0\\0&-m_{32}&1 \end{matrix}\right]U= \left[\begin{matrix} 1&0&0\\-m_{21}&1&0\\-m_{31}&-m_{32}&1 \end{matrix}\right]U \]
至此，我們已經成功對\(A\)進行了LU分解
\[L=\left[\begin{matrix} 1&0&0\\-m_{21}&1&0\\-m_{31}&-m_{32}&1 \end{matrix}\right], U=\left[\begin{matrix} a_{11}&a_{12}&a_{13}\\&a_{22}^{(1)}&a_{23}^{(1)}\\&&a_{33}^{(2)} \end{matrix}\right] \]
上述分析對於一般的\(n\)階矩陣都是適用的，但需要注意的是，LU分解是通過高斯消元法推導出來的，因此只有當\(A\)滿足順序主子式均大於零的時候，才能夠進行LU分解。上述過程涉及到一個小知識點：兩個下三角陣的乘積一定還是下三角陣。
在編程實現LU分解時，通過上述推導可以發現，在構造上三角陣\(U\)的時候，也可以同時得到下三角陣\(L\)。也就是說，當我們執行第\(k\)步高斯消元法，將\(a_{kk}^{(k-1)}\)下方的元素變成0時，需要先計算\(m_{ik}(i=k+1,\cdots,n)\)，此時就可以同時得到\(L\)的第\(k\)列元素\(-m_{ik}\)，而不需要單獨去求解\(L\)。
4.3.2 追趕法
在工程計算中，經常會遇到如下特殊的三對角線方程組：
\[Ax=\left[\begin{matrix} b_1&c_1&&&\\a_2&b_2&c_2&&\\ &\ddots&\ddots&\ddots&\\ &&a_{n-1}&b_{n-1}&c_{n-1}\\ &&&a_n&b_n \end{matrix}\right] \left[\begin{matrix} x_1\\x_2\\\vdots\\x_{n-1}\\x_n \end{matrix}\right]= \left[\begin{matrix} d_1\\d_2\\\vdots\\d_{n-1}\\d_n \end{matrix}\right]=d \]
所謂追趕法就是將LU分解應用到上述系數矩陣，從而得到的一個特殊解法。
首先對系數矩陣進行一個簡單的分析，如果對\(A\)進行高斯消元，可以發現所有\(c_i(1\le i\le n-1)\)的上方都是0，這意味着消元過程中，主對角線上方的三角形區域會完全保留在\(U\)中。
其次，主對角線下方只有\(a_i(2\le i\le n)\)需要消去，因此，第\(k\)步消元只需要計算\(m_{(k+1)k}\)，這意味着\(L\)下方只有緊挨着主對角線的元素不為0，其余均為0。
綜上分析，可以得到系數矩陣\(A\)的LU分解形式：
\[\left[\begin{matrix} b_1&c_1&&&\\a_2&b_2&c_2&&\\ &\ddots&\ddots&\ddots&\\ &&a_{n-1}&b_{n-1}&c_{n-1}\\ &&&a_n&b_n \end{matrix}\right]= \left[\begin{matrix} 1&&&&\\p_2&1&&&\\ &\ddots&\ddots&&\\ &&p_{n-1}&1&\\ &&&p_n&1 \end{matrix}\right] \left[\begin{matrix} q_1&c_1&&&\\&q_2&c_2&&\\ &&\ddots&\ddots&\\ &&&q_{n-1}&c_{n-1}\\ &&&&q_n \end{matrix}\right] \]
分析\(A\)的第\(k\)行如何通過右側兩個矩陣的乘積得到，根據矩陣分塊的知識，可以知道
\[\left[\begin{matrix} \cdots&\cdots&a_{k}&b_k&c_k&\cdots \end{matrix}\right]= \left[\begin{matrix} \cdots&\cdots&p_kq_{k-1}&p_kc_{k-1}+q_k&c_k&\cdots \end{matrix}\right] \]
上述省略號代表0，於是
\[a_k=p_kq_{k-1}\\ b_k=p_kc_{k-1}+q_k \]
可以驗證對於\(2\le k\le n\)，上述式子都是成立的，並且\(q_1=b_1\)。綜上可得追趕法計算公式：
\[\left\{\begin{aligned} q_1&=b_1\\ p_k&=\frac{a_k}{q_{k-1}}\\ q_k&=b_k-p_kc_{k-1} \end{aligned}\right.k=2,3,\dots,n \]
追趕法名稱的由來可能就是因為需要交替計算\(p_k,q_k\)，這相當於是\(p,q\)兩個下標之間在相互追趕。
求解\(Ly=d\)
\[\left[\begin{matrix} 1&&&&\\p_2&1&&&\\ &\ddots&\ddots&&\\ &&p_{n-1}&1&\\ &&&p_n&1 \end{matrix}\right]\left[\begin{matrix} y_1\\y_2\\\vdots\\y_{n-1}\\y_n \end{matrix}\right]= \left[\begin{matrix} d_1\\d_2\\\vdots\\d_{n-1}\\d_n \end{matrix}\right]\Rightarrow \left\{\begin{aligned} y_1&=d_1\\ p_2y_1+y_2&=d_2\\ \vdots&\\ p_ky_{k-1}+y_k&=d_k\\ \vdots&\\ p_ny_{n-1}+y_n&=d_n \end{aligned}\right. \]
因此
\[\left\{\begin{aligned} y_1&=d_1\\ y_k&=d_k-p_ky_{k-1} \end{aligned}\right.k=2,3,\cdots,n \]
求解\(Ux=y\)
\[\left[\begin{matrix} q_1&c_1&&&\\&q_2&c_2&&\\ &&\ddots&\ddots&\\ &&&q_{n-1}&c_{n-1}\\ &&&&q_n \end{matrix}\right] \left[\begin{matrix} x_1\\x_2\\\vdots\\x_{n-1}\\x_n \end{matrix}\right] = \left[\begin{matrix} y_1\\y_2\\\vdots\\y_{n-1}\\y_n \end{matrix}\right] \Rightarrow \left\{\begin{aligned} q_1x_1+c_1x_2&=y_1\\ \vdots&\\ q_kx_k+c_kx_{k+1}&=y_k\\ \vdots&\\ q_nx_n&=y_n \end{aligned}\right. \]
因此
\[\left\{\begin{aligned} x_n&=y_n/q_n\\ x_k&=(y_k-c_kx_{k+1})/q_k \end{aligned}\right.k=n-1,n-2,\cdots,1 \]
綜上可知，使用追趕法求解三對角線方程時，不需要真的去求解系數矩陣的LU分解，只需要根據上述的遞推公式就可以直接求問題的解。
