跳转至

高斯消元

引入

高斯消元法(Gauss–Jordan elimination)是求解線性方程組的經典算法,它在當代數學中有着重要的地位和價值,是線性代數課程教學的重要組成部分。

高斯消元法除了用於線性方程組求解外,還可以用於行列式計算、求矩陣的逆,以及其他計算機和工程方面。

消元法及高斯消元法思想

定義

消元法是將方程組中的一方程的未知數用含有另一未知數的代數式表示,並將其帶入到另一方程中,這就消去了一未知數,得到一解;或將方程組中的一方程倍乘某個常數加到另外一方程中去,也可達到消去一未知數的目的。消元法主要用於二元一次方程組的求解。

解釋

例一:利用消元法求解二元一次線性方程組:

\[ \begin{cases} 4x+y&=100 \\ x-y&=100 \end{cases} \]

解:將方程組中兩方程相加,消元 \(y\) 可得:

\(5x = 200\)

解得:

\(x = 40\)

\(x = 40\) 代入方程組中第二個方程可得:

\(y = -60\)

消元法理論的核心

消元法理論的核心主要如下:

  • 兩方程互換,解不變;

  • 一方程乘以非零數 \(k\),解不變;

  • 一方程乘以數 \(k\) 加上另一方程,解不變。

高斯消元法思想概念

德國數學家高斯對消元法進行了思考分析,得出瞭如下結論:

  • 在消元法中,參與計算和發生改變的是方程中各變量的係數;

  • 各變量並未參與計算,且沒有發生改變;

  • 可以利用係數的位置表示變量,從而省略變量;

  • 在計算中將變量簡化省略,方程的解不變。

高斯在這些結論的基礎上,提出了高斯消元法,首先將方程的增廣矩陣利用行初等變換化為行最簡形,然後以線性無關為準則對自由未知量賦值,最後列出表達方程組通解。

高斯消元五步驟法

解釋

高斯消元法在將增廣矩陣化為最簡形後對於自由未知量的賦值,需要掌握線性相關知識,且賦值存在人工經驗的因素,使得在學習過程中有一定的困難,將高斯消元法劃分為五步驟,從而提出五步驟法,內容如下:

  1. 增廣矩陣行初等行變換為行最簡形;

  2. 還原線性方程組;

  3. 求解第一個變量;

  4. 補充自由未知量;

  5. 列表示方程組通解。

利用實例進一步説明該算法的運作情況。

過程

例二:利用高斯消元法五步驟法求解線性方程組:

\[ \begin{cases} 2x_1+5x_3+6x_4&=9 \\ x_3+x_4&=-4 \\ 2x_3+2x_4&=-8 \end{cases} \]

增廣矩陣行(初等)變換為行最簡形

所謂增廣矩陣,即為方程組係數矩陣 \(A\) 與常數列 \(b\) 的並生成的新矩陣,即 \((A | b)\),增廣矩陣行初等變換化為行最簡形,即是利用了高斯消元法的思想理念,省略了變量而用變量的係數位置表示變量,增廣矩陣中用豎線隔開了係數矩陣和常數列,代表了等於符號。

\[ \left(\begin{matrix} 2 & 0 & 5 & 6 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 2 & 2 \end{matrix} \middle| \begin{matrix} 9 \\ -4 \\ -8 \end{matrix} \right) \]
\[ \xrightarrow{r_3-2r_2} \left(\begin{matrix} 2 & 0 & 5 & 6 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} \middle| \begin{matrix} 9 \\ -4 \\ 0 \end{matrix} \right) \]

化為行階梯形

\[ \xrightarrow{\frac{r_1}{2}} \left(\begin{matrix} 1 & 0 & 2.5 & 3 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} \middle| \begin{matrix} 4.5 \\ -4 \\ 0 \end{matrix} \right) \]
\[ \xrightarrow{r_1-r_2 \times 2.5} \left(\begin{matrix} 1 & 0 & 0 & 0.5 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{matrix} \middle| \begin{matrix} 14.5 \\ -4 \\ 0 \end{matrix} \right) \]

化為最簡形

還原線性方程組

\[ \begin{cases} x_1+0.5x_4 &= 14.5\\ x_3+x_4 &= -4 \\ \end{cases} \]

解釋

所謂的還原線性方程組,即是在行最簡形的基礎上,將之重新書寫為線性方程組的形式,即將行最簡形中各位置的係數重新賦予變量,中間的豎線還原為等號。

求解第一個變量

\[ \begin{cases} x_1 = -0.5x_4+14.5\notag \\ x_3 = -x_4-4\notag \end{cases} \]

解釋

即是對於所還原的線性方程組而言,將方程組中每個方程的第一個變量,用其他量表達出來。如方程組兩方程中的第一個變量 \(x_1\)\(x_3\)

補充自由未知量

\[ \begin{cases} x_1 = -0.5x_4+14.5 \\ x_2 = x_2 \\ x_3 = -x_4-4 \\ x_4 = x_4 \end{cases} \]

解釋

第 3 步中,求解出變量 \(x_1\)\(x_3\),從而説明了方程剩餘的變量 \(x_2\)\(x_4\) 不受方程組的約束,是自由未知量,可以取任意值,所以需要在第 3 步驟解得基礎上進行解得補充,補充的方法為 \(x_2 = x_2,x_4 = x_4\),這種解得補充方式符合自由未知量定義,並易於理解,因為是自由未知量而不受約束,所以只能自己等於自己。

列表示方程組的通解

\[ \begin{aligned} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix} &= \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix} x_2+ \begin{pmatrix} -0.5 \\ 0 \\ -1 \\ 1 \end{pmatrix} x_4 + \begin{pmatrix} 14.5 \\ 0 \\ -4 \\ 0 \end{pmatrix} \\ &= \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix} C_1+ \begin{pmatrix} -0.5 \\ 0 \\ -1 \\ 1 \end{pmatrix} C_2 + \begin{pmatrix} 14.5 \\ 0 \\ -4 \\ 0 \end{pmatrix} \end{aligned} \]

其中 \(C_1\)\(C_2\) 為任意常數。

解釋

即在第 4 步的基礎上,將解表達為列向量組合的表示形式,同時由於 \(x_2\)\(x_4\) 是自由未知量,可以取任意值,所以在解得右邊,令二者分別為任意常數 \(C_1\)\(C_2\),即實現了對方程組的求解。


行列式計算

解釋

\(N \times N\) 方陣行列式(Determinant)可以理解為所有列向量所夾的幾何體的有向體積

例如:

\[ \begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix} = 1 \]
\[ \begin{vmatrix} 1 & 2 \\ 2 & 1 \end{vmatrix} = -3 \]

行列式有公式

\[ \operatorname{det}(A)=\sum_{\sigma \in S_{n}} \operatorname{sgn}(\sigma) \prod_{i=1}^{n} a_{i, \sigma(i)} \]

其中 \(S_n\) 是指長度為 \(n\) 的全排列的集合,\(\sigma\) 就是一個全排列,如果 \(\sigma\) 的逆序對對數為偶數,則 \(\operatorname{sgn}(\sigma)=1\),否則 \(\operatorname{sgn}(\sigma)=−1\)

通過體積概念理解行列式不變性是一個非常簡單的辦法:

  • 矩陣轉置,行列式不變;

  • 矩陣行(列)交換,行列式取反;

  • 矩陣行(列)相加或相減,行列式不變;

  • 矩陣行(列)所有元素同時乘以數 \(k\),行列式等比例變大。

由此,對矩陣應用高斯消元之後,我們可以得到一個對角線矩陣,此矩陣的行列式由對角線元素之積所決定。其符號可由交換行的數量來確定(如果為奇數,則行列式的符號應顛倒)。因此,我們可以在 \(O(n^3)\) 的複雜度下使用高斯算法計算矩陣。

注意,如果在某個時候,我們在當前列中找不到非零單元,則算法應停止並返回 0。


實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
const double EPS = 1E-9;
int n;
vector<vector<double> > a(n, vector<double>(n));

double det = 1;
for (int i = 0; i < n; ++i) {
  int k = i;
  for (int j = i + 1; j < n; ++j)
    if (abs(a[j][i]) > abs(a[k][i])) k = j;
  if (abs(a[k][i]) < EPS) {
    det = 0;
    break;
  }
  swap(a[i], a[k]);
  if (i != k) det = -det;
  det *= a[i][i];
  for (int j = i + 1; j < n; ++j) a[i][j] /= a[i][i];
  for (int j = 0; j < n; ++j)
    if (j != i && abs(a[j][i]) > EPS)
      for (int k = i + 1; k < n; ++k) a[j][k] -= a[i][k] * a[j][i];
}

cout << det;

矩陣求逆

對於方陣 \(A\),若存在方陣 \(A^{-1}\),使得 \(A \times A^{-1} = A^{-1} \times A = I\),則稱矩陣 \(A\) 可逆,\(A^{-1}\) 被稱為它的逆矩陣。

給出 \(n\) 階方陣 \(A\),求解其逆矩陣的方法如下:

  1. 構造 \(n \times 2n\) 的矩陣 \((A, I_n)\)
  2. 用高斯消元法將其化簡為最簡形 \((I_n, A^{-1})\),即可得到 \(A\) 的逆矩陣 \(A^{-1}\)。如果最終最簡形的左半部分不是單位矩陣 \(I_n\),則矩陣 \(A\) 不可逆。

該方法的正確性證明需要用到較多線性代數的知識,限於篇幅這裏不再給出。感興趣的讀者可以自行查閲相關資料。

高斯消元法解異或方程組

異或方程組是指形如

\[ \begin{cases} a_{1,1}x_1 \oplus a_{1,2}x_2 \oplus \cdots \oplus a_{1,n}x_n &= b_1\\ a_{2,1}x_1 \oplus a_{2,2}x_2 \oplus \cdots \oplus a_{2,n}x_n &= b_2\\ \cdots &\cdots \\ a_{m,1}x_1 \oplus a_{m,2}x_2 \oplus \cdots \oplus a_{m,n}x_n &= b_1 \end{cases} \]

的方程組,其中 \(\oplus\) 表示「按位異或」(即 xor 或 C++ 中的 ^),且式中所有係數/常數(即 \(a_{i,j}\)\(b_i\))均為 \(0\)\(1\)

由於「異或」符合交換律與結合律,故可以按照高斯消元法逐步消元求解。值得注意的是,我們在消元的時候應使用「異或消元」而非「加減消元」,且不需要進行乘除改變係數(因為係數均為 \(0\)\(1\))。

注意到異或方程組的增廣矩陣是 \(01\) 矩陣(矩陣中僅含有 \(0\)\(1\)),所以我們可以使用 C++ 中的 std::bitset 進行優化,將時間複雜度降為 \(O(\dfrac{n^2m}{\omega})\),其中 \(n\) 為元的個數,\(m\) 為方程條數,\(\omega\) 一般為 \(32\)(與機器有關)。

參考實現:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
std::bitset<1010> matrix[2010];  // matrix[1~n]:增廣矩陣,0 位置為常數

std::vector<bool> GaussElimination(
    int n, int m)  // n 為未知數個數,m 為方程個數,返回方程組的解
                   // (多解 / 無解返回一個空的 vector)
{
  for (int i = 1; i <= n; i++) {
    int cur = i;
    while (cur <= m && !matrix[cur].test(i)) cur++;
    if (cur > m) return std::vector<bool>(0);
    if (cur != i) swap(matrix[cur], matrix[i]);
    for (int j = 1; j <= m; j++)
      if (i != j && matrix[j].test(i)) matrix[j] ^= matrix[i];
  }
  std::vector<bool> ans(n + 1, 0);
  for (int i = 1; i <= n; i++) ans[i] = matrix[i].test(0);
  return ans;
}

練習題