用消去法解線性方程

線性聯立方程組

  1. System of linear equations (線性聯立方程組): 雞兔同籠問題的延伸.
  2. coefficient 係數, variable 變數, constant term 常數項, solution 解.
  3. Q: 變數個數與方程式個數之間有什麼關係? 直覺猜測: 方程式太多 (稱為 overdetermined system), 可能會無解; 方程式太少 (稱為 underdetermined system), 可能會有很多組解.
  4. 結論: 不論方程式個數 (比起變數個數) 多寡, 都有可能無解, 也有可能無窮多組解 (Q: 用圖形舉例).
  5. 正確的直覺想法: 代數觀點 <--> 幾何意義

解線性聯立方程組的幾何解釋

解方程組所對應到的幾何問題:

代數觀點 幾何觀點
解線性聯立方程組 求超平面的交集
變數個數 空間維度
方程式個數 超平面個數
有幾組解? 交集有幾個點?

解集合所對應到的幾何圖形:

名稱 代數狀況 幾何狀況
inconsistent 無解 無交集
consistent 一組解 交集為一點
無限多組解 交集為 1-flat 或 2-flat 或 ...

方程式基本運算的幾何解釋

代數運算 幾何意義
兩方程式對調 看不出任何差別
方程式乘常數倍 原超平面
把另一方程式乘常數倍加入此方程式 過原來兩個超平面交集的另一超平面

Gaussian Elimination with Back Substitution

面對一組待解的線性方程組, 我們可以對其中的一個或一對方程式做上述三個簡單的動作, 不會影響 (改變) 答案. 這三個動作叫做 elementary row operations. 小時候解雞兔同籠的問題時, 所用的「加減消去法」, 所用的其實就是這三個運算.

為了少打 (寫) 一些字, 我們把方程組的係數及常數項拉出來 寫成一個矩陣, 稱為 augmented matrix. (反正變數名稱對計算過程沒有影響, 何苦每次都重抄一遍?) 例如要解:

          w + 3 x - 5 y + 2 z =  94
        4 w - 3 x +   y + 5 z =  58
        6 w - 2 x + 2 y + 4 z =  72
              2 x + 3 y - 7 z = -69

我們只需要看它的 augmented matrix

     /  1     3    -5     2   |  94 \
     |  4    -3     1     5   |  58 |
     |  6    -2     2     4   |  72 |
     \  0     2     3    -7   | -69 /

並且時時記得: 等一下對矩陣所做的運算, 其實只是在對原方程組做加減消去法, 沒什麼高深的學問.

做加減消去的目的何在? 希望把多一些係數化為 0, 讓方程組變得很簡單. 例如我們可以用第二列第三行的 1 把其他各列的第三行元素都 "做掉":

     / 21   -12     0    27   |  384 \ 
     |  4    -3     1     5   |   58 | 
     | -2     4     0    -6   |  -44 | 
     \-12    11     0   -22   | -243 /

這個動作稱為 選取第二列第三行元素做為 pivot. 第二列稱為 pivot row; 第三行稱為 pivot column. 當然如果選到的 pivot 元素不是 1, 不妨先將 pivot row 整列除以一個常數 (其實就是 pivot element 的值), 總之就是要先把 pivot 元素化為 1, 再來做加減消去. 例如從原來的矩陣出發, 如果當初我們不是選第二列第三行的 1 做為 pivot, 而是選第三列第二行的 -2, 那麼被 "做掉" 的, 就是其他列的第二行元素, 而矩陣就變成了:

     / 10     0    -2     8   |  202 \
     | -5     0    -2    -1   |  -50 |
     | -3     1    -1    -2   |  -36 |
     \  6     0     5    -3   |    3 /

作業: 請參考 「矩陣計算機 octave 與 rlab」 進入 octave 或 rlab 手動練習做 pivot 的動作.

我們希望重複上面 「選取 pivot 元素」 的動作很多次, 把矩陣的大部分元素都 "做掉". (其實只要能夠做掉大約一半不到的元素, 接下來的問題就很好解了.) 上面的例子目的只是在解釋何謂 pivoting; 但若要有規律地把矩陣的大部分元素都 "做掉", 下面的順序比較有規律, 有效率:

  1. 以第一列的頭 (第一行元素) 為 pivot, 消掉下面三列的第一行元素
  2. 以第二列的頭 (第二行元素) 為 pivot, 消掉下面兩列的第二行元素
  3. 以第三列的頭 (第三行元素) 為 pivot, 消掉下面一列的第三行元素
  4. 以第四列的頭 (第四行元素) 為 pivot, 只把自己變成 1 就好.

重點在於消去時不必太勤勞, 只要往下消就好, 不需要往上消. 最後矩陣的左下半全部變成 0, 對角線全部變成 1:

     /  1     3    -5     2   |   94 \
     |  0     1  -1.4   0.2   | 21.2 |
     |  0     0     1    -1   |  -17 |
     \  0     0     0     1   |    8 /

化為 row-echelon form 的過程當中, 如果遇到某列的頭為 0 怎麼辦? 先把下面 "頭不為 0" 的列對調上來再做 pivot. 不要去動上面已處理完的列, 否則就前功盡棄了.

不要忘記這個矩陣其實是方程組的係數, 如果把最開始被我們省略掉的變數名稱寫回來, 就知道原始的方程組已被我們用加減消去法改寫成:

     w + 3 x -   5 y +   2 z =   94
           x - 1.4 y + 0.2 z = 21.2
                   y -     z =  -17
                           z =    8

從最後一式馬上可以讀出 z = 8, 代回第三式可解得 y = -9, 再代回第二式解得 x = 7, 最後代回第一式解得 w = 12. 對一個矩陣重複做 elementary row operations, 消到只剩下右上角的這個演算法稱為 Gaussian elimination; 得到的矩陣稱為原矩陣的 row echelon form; 最後逐一往回代的動作稱為 back substitution.

作業 1: 請參考 「矩陣計算機 octave 與 rlab」 當中介紹的 lademo.m 或 lademo.r 程式, 自己呼叫 pivot 函數很多次, 一步步完成上題. 又, 用亂數產生矩陣, 給自己多出幾題. (例如 A = rand(5,5); b = rand(5,1); T = [A,b]; )

作業 2: 用 lema 驗證你的計算。 例如: perl lema 0.txt 你的系統內必須安裝有 perl 及 perl-tk。

Gauss-Jordan Elimination

做完 Gaussian elimination 但還沒有做 back substitution 時, 其實可以進一步由後往前再做 pivot 的動作:

  1. 以第四列第四行元素為 pivot 消去上面三列的最右邊那個 (第四行) 元素, 得到: 
         /  1     3    -5     0   |   78 \
     |  0     1  -1.4     0   | 19.6 |
     |  0     0     1     0   |   -9 |
     \  0     0    -0     1   |    8 /
  2. 以第三列第三行元素為 pivot 消去上面兩列的第三行元素.
  3. 以第二列第二行元素為 pivot 消去上面一列的第二行元素.

最後矩陣左邊變成單位矩陣, 而右邊就是答案, 連 back substitution 都不必做. 這個演算法稱為 Gauss-Jordan elemination; 算出來的矩陣稱為原矩陣的 reduced row echelon form.

     w                       =   12
           x                 =    7
                   y         =   -9
                           z =    8

Singularity/Degeneracy

有很多數學公式都有 "極少發生/極不可能發生" 的例外狀況, 這種狀況稱為 singularity 或 degeneracy (形容詞為 singular/degenerate). 例:

  1. 若 x 為未知數, a, b 為已知數, 則 a * x = b 的解為,x = b/a。 但是如果 a=0 呢?
  2. 不管 x 是多少, x^0 = 1。 但是如果 x=0 呢?
  3. y = a*x^2 + b*x + c 的函數圖形是一條拋物線。 但是如果 a=0 呢?

以解聯立方程組而言, 如果你的方程組係數是用亂數產生的, 那麼上兩節的方法 "幾乎" 可以解決所有問題; 但是如果遇到鬼, 或是中了樂透頭獎 (意思是一樣的 -- 機率近乎 0), 也會碰到 singularity: 記得我們選的 pivot 元素不可為 0, 如果是 0 的話, 要到下面調一列頭不為 0 的列上來; 那如果下面所有列的頭都是 0 呢? 那就是遇到鬼了, 等一下會看到本題有無窮多組解或無解, 就是不會像平常一樣恰有一組解. 此時我們往右退一格, 總之要選 「剩下所有方程式當中能夠調到最左上角的非 0 元素」 當做 pivot.

只要在解題過程中遇到一次 我們就說這題是 singular/degenerate. 當然也有極極小的可能遇到第二次, ... 那就是 singularity 當中的 singularity, ... Singular 的問題, 如果最下面幾式都是 0 = 0 則有無窮多組解, 表示超平面的交集不只一個點, 是一個 dimension 大於或等於 1 的 flat. 從 row-echelon form 開始往回解的過程當中, 必會遇到變數太多, 條件不足的狀況, 此時用參數式表達. Singular 的問題, 如果最下面幾式出現 0 = 2.7 之類不合理的狀況, 表示超平面沒有交集, 此題無解. 例:

         - x  + 2 y             =  3
         2 x  - 4 y +   z + 3 w = -4
           x  - 2 y + 3 z + 9 w =  3
        -2 x  + 4 y + 2 z + 6 w = 10

將 augmented matrix 化為 reduced row echelon form 之後, 得到:

        / 1 -2  0  0 | -3 \
        | 0  0  1  3 |  2 |
        | 0  0  0  0 |  0 |
        \ 0  0  0  0 |  0 /

作業: 用 rndle 隨機產生一些題目, 再用 lema 解題。 例如:

        ./rndle -m 6 -n 6 -r 4 > test.le
        ./lema -f 12x24 test.le

產生一個大小為 6x6, 秩為 4 的題目, 並用 "12x24" 的字形顯示。 要看有那些不同大小的字形可用, 可以下 xlsfonts | grep '[0-9]x[0-9]'

本節適用於 lema 的練習題: 1.le, 2.le, 3.le, 4.le, 5.le

將 Gauss-Jordan Elimination 應用於矩陣

用 Gauss-Jordan Elimination 可以求一個 n*n 方陣 A 的反矩陣 (如果它有反矩陣的話): 要求 B 使得 A * B = I 其實可以看成是 n 題 線性聯立方程組的問題. 更好的是, 這 n 題的運算過程當中, 所用到的 elementary row operation (基本列運算) 一模一樣, 所以可以乾脆同時做: 把 A 與 I 並排, 用 Gauss-Jordan 把左邊的 A 變成 I 的同時, 右邊的 I 會自動變成 A 的反矩陣.

Gauss-Jordan elemination 不一定要用在 augmented matrix 上; 也就是說, 它不一定要與解聯立方程組有關. 如果把它用在普通的矩陣 A 上, 所得到的矩陣還是叫做 A 的 reduced row echelon form. 扣除掉此矩陣當中全為零的那幾列, 剩下來的所有列向量 (即所有「非零向量」的列向量) 構成 A 的 row space 的一組基底. 這些非零向量的個數, 稱為 rank of A A 的秩.

existence of inverse, singularity, existence of solutions, rank