高斯消去法

維基百科,自由的百科全書
跳轉到: 導覽搜尋
Confusion grey.svg
提示:本條目的主題不是高斯-若爾當消元法

數學上,高斯消去法(或譯:高斯消元法)(英語:Gaussian Elimination),是線性代數中的一個演算法,可用來為線性方程組求解,求出矩陣的秩,以及求出可逆方陣逆矩陣。當用於一個矩陣時,高斯消去法會產生出一個「行梯陣式」。

歷史[編輯]

該方法以數學家卡爾·高斯命名,但最早出現於中國古籍《九章算術》,成書於約公元前150年。[1]

例子[編輯]

高斯消去法可用來找出下列方程組的解或其解的限制:

2x + y - z = 8 \quad (L_1)
-3x - y + 2z = -11 \quad (L_2)
-2x + y + 2z = -3 \quad (L_3)

這個演算法的原理是:

首先,要將L_1以下的等式中的x消除,然後再將L_2以下的等式中的y消除。這樣可使整個方程組變成一個三角形似的格式。之後再將已得出的答案一個個地代入已被簡化的等式中的未知數中,就可求出其餘的答案了。

在剛才的例子中,我們將\begin{matrix}\frac{3}{2}\end{matrix} L_1L_2相加,就可以將L_2中的x消除了。然後再將L_1L_3相加,就可以將L_3中的x消除。

我們可以這樣寫:

L_2 + \frac{3}{2}L_1 \rightarrow L_2
L_3 + L_1 \rightarrow L_3

結果就是:

2x + y - z = 8 \,
\frac{1}{2}y + \frac{1}{2}z = 1 \,
2y + z = 5 \,

現在將-4L_2L_3相加,就可將L_3中的y消除:

L_3 + -4L_2 \rightarrow L_3

其結果是:

2x + y - z = 8 \,
\frac{1}{2}y + \frac{1}{2}z = 1 \,
-z = 1 \,

這樣就完成了整個演算法的初步,一個三角形的格式(指:變數的格式而言,上例中的變數各為3,2,1個)出現了。

第二步,就是由尾至頭地將已知的答案代入其他等式中的未知數。第一個答案就是:

z = -1 \quad (L_3)

然後就可以將z代入L_2中,立即就可得出第二個答案:

y = 3 \quad (L_2)

之後,將zy代入L_1之中,最後一個答案就出來了:

x = 2 \quad (L_1)

就是這樣,這個方程組就被高斯消去法解決了。

這種演算法可以用來解決所有線性方程組。即使一個方程組不能被化為一個三角形的格式,高斯消去法仍可找出它的解。例如,如果在第一步化簡後,L_2L_3中沒有出現任何y,沒有三角形的格式,照著高斯消去法而產生的格式仍是一個行梯陣式。這情況之下,這個方程組會有超過一個解,當中會有至少一個變數作為答案。每當變數被鎖定,就會出現一個解。

通常人或電腦在應用高斯消去法的時候,不會直接寫出方程組的等式來消去未知數,反而會使用矩陣來計算。以下就是使用矩陣來計算的例子:


\begin{pmatrix}
2 & 1 & -1 & 8 \\
-3 & -1 & 2 & -11 \\
-2 & 1 & 2 & -3
\end{pmatrix}

跟著以上的方法來運算,這個矩陣可以轉變為以下的樣子:


\begin{pmatrix}
2 & 1 & -1 & 8 \\
0 & \frac{1}{2} & \frac{1}{2} & 1 \\
0 & 0 & -1 & 1
\end{pmatrix}

這矩陣叫做「行梯陣式」。

最後,可以利用同樣的演算法產生以下的矩陣,便可把所得出的解或其限制簡明地表示出來:


\begin{pmatrix}
1 & 0 & 0 & 2 \\
0 & 1 & 0 & 3 \\
0 & 0 & 1 & -1
\end{pmatrix}

最後這矩陣叫做「簡化行梯陣式」,亦是高斯-約當消去法指定的步驟。

其他應用[編輯]

找出逆矩陣[編輯]

高斯消去法可以用來找出一個可逆矩陣逆矩陣。設A為一個n \times n的矩陣,其逆矩陣可被兩個分塊矩陣表示出來。將一個n \times n 單位矩陣放在A的右手邊,形成一個n \times 2n的分塊矩陣B = [A, I]。經過高斯消去法的計算程序後,矩陣B的左手邊會變成一個單位矩陣I,而逆矩陣A^{-1}會出現在B的右手邊。

假如高斯消去法不能將A化為三角形的格式,那就代表A是一個不可逆的矩陣。

應用上,高斯消去法極少被用來求出逆矩陣。高斯消去法通常只為線性方程組求解[2]

計出秩的基本演算法[編輯]

高斯消去法可應用在任何m \times n矩陣A。在不可減去某數的情況下,我們都只有跳到下一行。以一個6 \times 9的矩陣作例,它可以變化為一個行梯陣式:


\begin{pmatrix}
1 & * & 0 & 0 & * & * & 0 & * & 0 \\
0 & 0 & 1 & 0 & * & * & 0 & * & 0 \\
0 & 0 & 0 & 1 & * & * & 0 & * & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & * & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{pmatrix}

而矩陣中的*'是一些數字。這個梯陣式的矩陣T會有一些關於A的資訊:

  • A是5,因為T有5列非0的列;
  • A的列的向量空間,可從A的第1、3、4、7和9列中得知,其數值在矩陣T之中;
  • 矩陣中的*'表示了A的列可怎樣寫為列中的數的組合。

分析[編輯]

高斯消去法演算法複雜度O(n3);這就是說,如果係數矩陣的是n × n,那麼高斯消去法所需要的計算量大約與n3成比例

高斯消去法可以用在電腦中來解決數千條等式及未知數。不過,如果有過百萬條等式時,這個演算法會十分費時。一些極大的方程組通常會用迭代法來解決。亦有一些方法特地用來解決一些有特別排列的係數的方程組。

高斯消去法可用在任何中。

高斯消去法對於一些矩陣來說是穩定的。對於普遍的矩陣來說,高斯消去法在應用上通常也是穩定的,不過亦有例外。[3]

偽代碼[編輯]

高斯消去法的其中一種偽代碼

i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:    // 从第i行开始,找出第j列中的最大值(i、j值应保持不变)  
  maxi := i
  for k := i+1 to m do
    if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi := k   // 使用交换法找出最大值(绝对值最大)
    end if
  end for
  if A[maxi,j] ≠ 0 then  // 判定找到的绝对值最大值是否为零:若不为零就进行以下操作;若为零则说明该列第(i+1)行以下(包括第(i+1)行)均为零,不需要再处理,直接跳转至第(j+1)列第(i+1)行
    swap rows i and maxi, but do not change the value of i   // 将第i行与找到的最大值所在行做交换,保持i值不变(i值记录了本次操作的起始行)
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry in row i by A[i,j]    // 将交换后的第i行归一化(第i行所有元素分别除以A[i,j])
    Now A[i,j] will have the value 1.
    for u := i+1 to m do    // 第j列中,第(i+1)行以下(包括第(i+1)行)所有元素都减去A[i,j]
      subtract A[u,j] * row i from row u
      Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end for
    i := i + 1   
  end if
  j := j + 1  // 第j列中,第(i+1)行以下(包括第(i+1)行)所有元素均为零。移至第(j+1)列,从第(i+1)行开始重复上述步骤。
end while

這個演算法和之前談到的有點不同,它由絕對值最大的部分開始做起,這樣可以改善演算法的穩定性。本演算法由左至右地計算,每作出以下三個步驟,才跳到下一列和下一行:

  1. 定出每列的絕對值最大的一個非0的數,將第一行的值與該行交換,使得第一行擁有這一列的最大值;
  2. 將第一列的數字除以該數,使得該行的第一個數成為1;
  3. 對每一行減去第一行乘以每一行的第一個數,使得每一行的第一個數變為0。

所有步驟完成後,這個矩陣會變成一個行梯矩陣,再用代入法就可以求解該方程組。


隨著多核處理器的日益普及,現在的程式設計師可以利用執行緒級並行高斯消元演算法來提高計算的速度。內存共享模式(而不是消息交換模式)的偽代碼如下所示:

// Note this code sample has been mangled (missing attr init/scope? bad gauss() indentation?)...
// What is original source?  Can we get valid C or else simplified pseudo code instead of this hybrid?
 
 void parallel(int num_threads,int matrix_dimension)
 {
  int i;
  for(i=0;i<num_threads;i++)
    create_thread(&threads[i],i);
  pthread_attr_destroy(&attr); // Free attribute and wait for the other threads
  for(i=0;i<p;i++)
    pthread_join(threads[i],NULL);
 }
 void *gauss(int thread_id)
 {
  int i,k,j;
  for(k=0;k<matrix_dimension-1;k++)
  {
    if(thread_id==(k%num_thread))  //interleaved-row work distribution
    {
      for(j=k+1;j<matrix_dimension;j++)
        M[k][j]=M[k][j]/M[k][k];
      M[k][k]=1;
    }
    barrier(num_thread,&mybarrier);      //wait for other thread finishing this round
    for(i=k+1;i<matrix_dimension;i=i+1)
        if(i%p==thread_id)
        {
           for(j=k+1;j<matrix_dimension;j++)
             M[i][j]=M[i][j]-M[i][k]*M[k][j];
           M[i][k]=0;
        }
    barrier(num_thread,&mybarrier);
  }
  return NULL;
 }
 void barrier(int num_thread, barrier_t * mybarrier)
 {
   pthread_mutex_lock(&(mybarrier->barrier_mutex));
   mybarrier->cur_count++;
   if(mybarrier->cur_count!=num_thread)
        pthread_cond_wait(&(mybarrier->barrier_cond),&(mybarrier->barrier_mutex));
   else
   {
       mybarrier->cur_count=0;
       pthread_cond_broadcast(&(mybarrier->barrier_cond));
   }
   pthread_mutex_unlock(&(mybarrier->barrier_mutex));
 }

參考文獻[編輯]

  1. ^ 第八章 方程//九章算術》. 150bc [2007年12月25日]. 
  2. ^ Atkinson, 1989年,第514頁
  3. ^ Golub and Van Loan, §3.4.6
  • Atkinson, Kendall A. An Introduction to Numerical Analysis, 第二版, John Wiley & Sons, New York, 1989年 ISBN 0-471-50023-2
  • Golub, Gene H., and Van Loan, Charles F. Matrix computations, 第三版, Johns Hopkins, Baltimore, 1996年 ISBN 0-8018-5414-8
  • Lipschutz, Seymour, and Lipson, Mark. Schaum's Outlines: Linear Algebra, Tata McGraw-hill edition.Delhi 2001年, 第69-80頁

參見[編輯]

外部連接[編輯]