高斯消去法

線性代數中求各種解的一個算法

高斯消去法(英語:Gaussian Elimination)是线性代数中的一个算法,可以把矩阵转化为行阶梯形矩阵[1]高斯消去法可用來為線性方程組求解,求出矩陣的秩,以及求出可逆方陣逆矩陣

历史

编辑

该方法以数学家卡尔·高斯命名,但最早出现于中国古籍《九章算术》,成书于约公元前150年。[2]

例子

编辑

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

這個算法的原理是:

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

在剛才的例子中,我們將  相加,就可以將 中的 消除了。然後再將  相加,就可以將 中的 消除。

我們可以這樣寫:

 
 

結果就是:

 
 
 

現在將  相加,就可將 中的 消除:

 

其結果是:

 
 
 

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

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

 

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

 

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

 

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

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

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

 

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

 

這矩陣叫做「行阶梯形矩阵」。

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

 

最後這矩陣叫做「簡化行阶梯形矩阵」,亦是高斯-若爾當消元法指定的步驟。[3]

其他應用

编辑

找出逆矩陣

编辑

高斯消去法可以用來找出一個可逆矩陣逆矩陣。設 為一個 的矩陣,其逆矩陣可被兩個分塊矩陣表示出來。將一個  單位矩陣放在 的右手邊,形成一個 的分塊矩陣 。經過高斯消去法的計算程序後,矩陣 的左手邊會變成一個單位矩陣 ,而逆矩陣 會出現在 的右手邊。

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

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

計出秩和基底的基本算法

编辑

高斯消去法可應用在任何 矩陣 。在不可減去某數的情況下,我們都只有跳到下一行。以一個 的矩陣作例,它可能可以變化為一個行阶梯形矩阵:

 

而矩陣中的*是一些數字。這個行阶梯形矩阵 會有一些關於 的資訊:

  •  是5,因為 有5行非0的行;
  •  列空間Col(A),將以 的第1、3、4、7和9行為基底,就是那些在矩陣 之中擁有主元的行;
  • 矩陣中的*表示了 的該行如何寫為上述基底的線性組合

分析

编辑

高斯消去法的算法复杂度O(n3);这就是说,如果系數矩阵的是n × n,那么高斯消去法所需要的计算量大约与n3成比例

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

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

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

伪代码

编辑

高斯消去法的其中一种伪代码

 i = 1
 j = 1
 while (i  m and j  n) do
   Find pivot in column j, starting in row i    // 从第i行(row)开始,找出第j列(column )中的最大值(i、j值应保持不变)  #台湾与大陆的列、行定义相反。台湾列为row行为column,大陆列为column行为row。
   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],直到第j行的i+1列以後元素均為零
       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。

所有步骤完成后,这个矩阵会变成一个行阶梯形矩阵,再用代入法就可以求解该方程组。


随着多核处理器的日益普及,现在的程序员可以利用线程级并行高斯消元算法来提高计算的速度。内存共享模式(而不是消息交换模式)的伪代码如下所示:

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][j]*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. ^ Steven J. Leon. Linear Algebra With Applications. 培生教育. 2014: 第13頁. ISBN 9781292025148. 
  2. ^ 第八章 方程. 《九章算术》. 150bc [2007年12月25日]. (原始内容存档于2009年4月19日). 
  3. ^ Steven J. Leon. Linear Algebra With Applications. 培生教育. 2014: 第17頁. ISBN 9781292025148. 
  4. ^ Atkinson, 1989年,第514頁
  5. ^ Golub and Van Loan, §3.4.6
  • Atkinson, Kendall A. An Introduction to Numerical Analysis, 第二版, John Wiley & Sons, New York, 1989年 ISBN 978-0-471-50023-0
  • Golub, Gene H., and Van Loan, Charles F. Matrix computations, 第三版, Johns Hopkins, Baltimore, 1996年 ISBN 978-0-8018-5414-9
  • Lipschutz, Seymour, and Lipson, Mark. Schaum's Outlines: Linear Algebra, Tata McGraw-hill edition.Delhi 2001年, 第69-80頁

參見

编辑

外部链接

编辑