高斯消元法(Gauss Elimination)【超詳解&模板】


高斯消元法,是線性代數中的一個算法,可用來求解線性方程組,並可以求出矩陣的秩,以及求出可逆方陣的逆矩陣。
高斯消元法的原理是:
若用初等行變換將增廣矩陣 化為 ,則AX = B與CX = D是同解方程組。

所以我們可以用初等行變換把增廣矩陣轉換為行階梯陣,然后回代求出方程的解。

 

1、線性方程組

                          image

1)構造增廣矩陣,即系數矩陣A增加上常數向量b(A|b)

                          image

2)通過以交換行、某行乘以非負常數和兩行相加這三種初等變化將原系統轉化為更簡單的三角形式(triangular form)

     注:這里的初等變化可以通過系數矩陣A乘上初等矩陣E來實現

                         image

3)從而得到簡化的三角方陣組,注意它更容易解

                         image

4)這時可以使用向后替換算法(Algorithm for Back Substitution)求解得

    z=-4/-4=1,  y=4-2z=4-2=2,  x= (1-y-z)/2=(1-2-1)/2=-1

 

總結上面過程,高斯消元法其實就是下面非常簡單的過程

                            原線性方程組       ——>       高斯消元法     ——> 下三角或上三角形式的線性方程組           ——>  前向替換算法求解(對於上三角形式,采用后向替換算法)

補充1:

 

高斯-若爾當消元法(Gauss-Jordan Elimination

 

相對於高斯消元法,高斯-若爾當消元法最后的得到線性方程組更容易求解,它得到的是簡化行列式。其轉化后的增高矩陣形式如下,因此它可以直接求出方程的解,而無需使用替換算法。但是,此算法的效率較低。

                             image

 

例子如下:

image          解為image

個人感覺區別就是對每行進行了歸一化處理

補充2:

介紹了最基本的高斯消元法,現在看看應用於實際問題的實用算法

因為實際應用中,我們總是利用計算機來分析線性系統,而計算機中以有限的數來近似無限的實數,因此產生舍入誤差(roundoff error),進而對解線性系統產生很多影響。

 

一個t位(即精度為t)以image為基的浮點數的表達形式為:imageimage。對於一個實數x,其浮點近似值image為最接近x的浮點數,必要時進行近似image

例1:對2位以10為基的浮點算法,image

例2:同樣考慮imageimage

 

以下面系統為例,看看在高斯消元中采用浮點算法會產生什么效果。

                                                                           image

當以精確解法時,通過將第一行乘以m=89/47,並從第二行中減去得到image,進而利用后向替換算法得x=1,y=-1。

當以3位以10為基的浮點算法時,乘子變為image,因為image,因此第一步高斯消元后得

image。此時,因為不能將第2行第1列位置變為0,所以不能將其三角化。從而,我們只能接受將這個位置值賦為0,而不管其實際浮點值。因此,3位浮點高斯消元的結果為image,后向算法計算結果為image

盡管無法消除近似誤差的影響,可以采用一些技術來盡量減小這類機器誤差。部分主元消元法在高斯消元的每一步,都選擇列上最大值為軸(通過行變換將其移動)。

下面給出列主元消去的代碼(所謂列主元消去法是在系數矩陣中按列選取元素絕對值得最大值作為主元素,然后交換所在行與主元素所在行的位置,再按順序消去法進行消元。)
 1 列選主元消元法
 2 /*
 3 *Gauss.cpp
 4 *功能:列選主元消元法 
 5 */
 6 #include<stdio.h>
 7 #include"Gass.h"
 8 
 9 //高斯消元法(列選主元)
10 void Gauss (double a[][MAXNUM],int n)
11 {
12     int i,j;
13     SelectColE(a,n);   //列選主元並消元成上三角
14     //回代求解
15     printf("上三角的結果\n");
16     printM(a,3);
17     for(i=n;i>=1;i--)
18     {
19         for(j=i+1;j<=n;j++)
20             a[i][n+1]-=a[i][j]*a[j][n+1];
21         a[i][n+1]/=a[i][i];
22     }
23     return ;
24 }
25 //選擇列主元並進行消元
26 void SelectColE(double a[][MAXNUM],int n)
27 {
28     int i,j,k,maxRowE;
29     double temp; //用於記錄消元時的因數
30     for(j=1;j<=n;j++)
31     {
32         maxRowE=j;
33         for(i=j;i<=n;i++)
34             if(fabs(a[i][j])>fabs(a[maxRowE][j]))
35                 maxRowE = i;
36         if(maxRowE!=j)
37             swapRow(a,j,maxRowE,n);   //與最大主元所在行交換
38         //消元
39         for(i=j+1;i<=n;i++)
40         {
41             temp =a[i][j]/a[j][j];
42             for(k=j;k<=n+1;k++)
43                 a[i][k]-=a[j][k]*temp;
44         }
45         printf("第%d列消元后:\n",j);
46         printM(a,3);
47     }
48     return;
49 }
50 void swapRow(double a[][MAXNUM],int m,int maxRowE,int n)
51 {
52     int k;
53     double temp;
54     for(k=m;k<=n+1;k++)
55     {
56         temp = a[m][k];
57         a[m][k] = a[maxRowE][k];
58         a[maxRowE][k] = temp;
59     }
60     return ;
61 }
62 
63 //測試函數
64 int main()
65 {
66     double a[4][MAXNUM];
67     
68     int i,n,j;
69 
70     a[1][1] = 0.001; a[1][2]=2.000;a[1][3]=3.000;a[1][4]=1.000;
71     a[2][1] = -1.000;a[2][2]=3.712;a[2][3]=4.623;a[2][4]=2.000;
72     a[3][1] = -2.000;a[3][2]=1.070;a[3][3]=5.643;a[3][4]=3.000;
73     Gauss (a,3);
74     for(i=1;i<=3;i++)
75         printf("X%d=%f\n",i,a[i][4]);
76     return 0;
77 }
78 //輸出矩陣
79 void printM(double a[][MAXNUM], int n)
80 {
81    int i,j;
82    for(i=1;i<=n;i++)
83     {
84         for(j=1;j<=n+1;j++)
85         {
86             printf("%f ,",a[i][j]);
87         }
88         printf("\n");
89     }
90 }
測試結果:

2、逆矩陣


 

 

下面來說說高斯消元法在編程中的應用。

首先,先介紹程序中高斯消元法的步驟:
(我們設方程組中方程的個數為equ,變元的個數為var,注意:一般情況下是n個方程,n個變元,但是有些題目就故意讓方程數與變元數不同)

1. 把方程組轉換成增廣矩陣。

2. 利用初等行變換來把增廣矩陣轉換成行階梯陣。
枚舉k從0到equ – 1,當前處理的列為col(初始為0) ,每次找第k行以下(包括第k行),col列中元素絕對值最大的列與第k行交換。如果col列中的元素全為0,那么則處理col + 1列,k不變。

3. 轉換為行階梯陣,判斷解的情況。

① 無解
當方程中出現(0, 0, …, 0, a)的形式,且a != 0時,說明是無解的。

② 唯一解
條件是k = equ,即行階梯陣形成了嚴格的上三角陣。利用回代逐一求出解集。

③ 無窮解。
條件是k < equ,即不能形成嚴格的上三角形,自由變元的個數即為equ – k,但有些題目要求判斷哪些變元是不缺定的。
    這里單獨介紹下這種解法:
首先,自由變元有var - k個,即不確定的變元至少有var - k個。我們先把所有的變元視為不確定的。在每個方程中判斷不確定變元的個數,如果大於1個,則該方程無法求解。如果只有1個變元,那么該變元即可求出,即為確定變元。

以上介紹的是求解整數線性方程組的求法,復雜度是O(n3)。浮點數線性方程組的求法類似,但是要在判斷是否為0時,加入EPS,以消除精度問題。

實例代碼:

有多組測試數據。每組測試數據先輸入一個整數n,表示方陣的階。然后下面輸入n階方陣。輸出其逆矩陣。若無逆矩陣,則輸出No inverse matrix。

 

  1 #include <iostream>
  2 #include <cmath>
  3 #include <algorithm>
  4 
  5 using namespace std;
  6 
  7 const double eps = 1e-6;
  8 
  9 bool is_zero( const double num )//用於判斷有無逆矩陣
 10 {
 11     return fabs(num) < eps;
 12 }
 13 
 14 void create( double ** & matrix, const int n )
 15 {
 16     matrix = new double* [n];
 17     for ( int i = 0; i < n; ++i )
 18         matrix[i] = new double[n];
 19 }
 20 
 21 void input ( double ** matrix, const int n )
 22 {
 23     for ( int i = 0; i < n; ++i )
 24     {
 25         for ( int j  = 0; j < n; ++ j )
 26             cin >> matrix[i][j];
 27     }
 28 }
 29 
 30 bool inverse ( double ** matrix1, double ** matrix2, const int n )
 31 {
 32     int i, j;
 33     for ( i = 0; i < n; ++ i )//初始化一個單位矩陣
 34     {
 35         for ( j = 0; j < n; ++ j )
 36         {
 37             if ( i == j )
 38                 matrix2[i][j] = 1;
 39             else
 40                 matrix2[i][j] = 0;
 41         }
 42     }
 43     for ( i = 0; i < n; ++i )
 44     {
 45         int rowmaxpos = i;
 46         for ( j = i + 1; j < n; ++j )
 47         {
 48             if ( matrix1[i][j] > matrix1[i][rowmaxpos] )
 49                 rowmaxpos = j;
 50         }
 51         for ( j = i; j < n; ++ j )//按從大到小的順序排列矩陣
 52         {
 53             swap( matrix1[j][rowmaxpos], matrix1[j][i]);
 54             swap( matrix2[j][rowmaxpos], matrix2[j][i]);
 55         }
 56         if ( !is_zero(matrix1[i][i]) )
 57         {
 58             int divisor = matrix1[i][i];
 59             for ( j = i; j < n; ++ j )//歸一化矩陣
 60             {
 61                 matrix1[i][j] /= divisor;
 62                 matrix2[i][j] /= divisor;
 63             }
 64             for ( j = i + 1; j < n; ++ j )//高斯消元法處理行列式
 65             {
 66                 int multiple = matrix1[j][i];
 67                 for ( int k = i; k < n; ++ k )
 68                 {
 69                     matrix1[i][j] -= matrix1[i][k] * multiple;
 70                     matrix2[i][j] -= matrix2[i][k] * multiple;
 71                 }
 72             }
 73         }
 74         else
 75             return false;
 76     }
 77     return true;
 78 }
 79 
 80 void output( double ** matrix, const int n )
 81 {
 82     for ( int i = 0; i < n; ++i )
 83     {
 84         for ( int j = 0; j < n; ++ j )
 85             cout << matrix[i][j] << ' ';
 86         cout<<endl;
 87     }
 88 }
 89 
 90 void destroy( double ** matrix, const int n )
 91 {
 92     for ( int i = 0; i < n; ++ i )
 93         delete [] matrix[i];
 94     delete [] matrix;
 95 }
 96 
 97 int main()
 98 {
 99     int n;
100     double ** matrix1;
101     double ** matrix2;
102     while ( cin >> n )
103     {
104         create( matrix1, n );
105         create( matrix2, n );
106         input( matrix1, n);
107         if ( inverse(matrix1, matrix2, n) )
108             output( matrix2, n );
109         else
110             cout << "No inverse matrix" << endl;
111         destroy( matrix1, n );
112         destroy( matrix2, n );
113     }
114     return 0;
115 }

 

擴展1:

 

利用矩陣的初等行變換也可以判斷一個矩陣是否可逆,即分塊矩陣(A︱E)經過初等行變換,原來A的位置不能變換為單位陣E,那么A不可逆。


擴展2:

 

利用矩陣初等行變換解矩陣方程

對於一般的矩陣方程,我們可以先求自變量系數矩陣的逆,然后乘以結果矩陣即可得到自變量矩陣

ps:最后來點福利:

 

下面是幾道OJ上的高斯消元法求解線性方程組的題目:

POJ 1222 EXTENDED LIGHTS OUT
http://acm.pku.edu.cn/JudgeOnline/problem?id=1222
POJ 1681 Painter's Problem
http://acm.pku.edu.cn/JudgeOnline/problem?id=1681
POJ 1753 Flip Game
http://acm.pku.edu.cn/JudgeOnline/problem?id=1753
POJ 1830 開關問題
http://acm.pku.edu.cn/JudgeOnline/problem?id=1830

POJ 3185 The Water Bowls

http://acm.pku.edu.cn/JudgeOnline/problem?id=3185
開關窗戶,開關燈問題,很典型的求解線性方程組的問題。方程數和變量數均為行數*列數,直接套模板求解即可。但是,當出現無窮解時,需要枚舉解的情況,因為無法判斷哪種解是題目要求最優的。

POJ 2947 Widget Factory
http://acm.pku.edu.cn/JudgeOnline/problem?id=2947
求解同余方程組問題。與一般求解線性方程組的問題類似,只要在求解過程中加入取余即可。
注意:當方程組唯一解時,求解過程中要保證解在[3, 9]之間。

POJ 1166 The Clocks
http://acm.pku.edu.cn/JudgeOnline/problem?id=1166
經典的BFS問題,有各種解法,也可以用逆矩陣進行矩陣相乘。
但是這道題用高斯消元法解決好像有些問題(困擾了我N天...持續困擾中...),由於周期4不是素數,故在求解過程中不能進行取余(因為取余可能導致解集變大),但最后求解集時,還是需要進行取余操作,那么就不能保證最后求出的解是正確的...在discuss里提問了好幾天也沒人回答...希望哪位路過的大牛指點下~~

POJ 2065 SETI
http://acm.pku.edu.cn/JudgeOnline/problem?id=2065
同樣是求解同余方程組問題,由於題目中的p是素數,可以直接在求解時取余,套用模板求解即可。(雖然AC的人很少,但它還是比較水的一道題,)

POJ 1487 Single-Player Games
http://acm.pku.edu.cn/JudgeOnline/problem?id=1487
很麻煩的一道題目...題目中的敘述貌似用到了編譯原理中的詞法定義(看了就給人不想做的感覺...)
解方程組的思想還是很好看出來了(前提是通讀題目不下5遍...),但如果把樹的字符串表達式轉換成方程組是個難點,我是用棧 + 遞歸的做法分解的。首先用棧的思想求出該結點的孩子數,然后遞歸分別求解各個孩子。
這題解方程組也與眾不同...首先是求解浮點數方程組,要注意精度問題,然后又詢問不確定的變元,按前面說的方法求解。
一頓折騰后,這題居然寫了6000+B...而且囧的是巨人C++ WA,G++ AC,可能還是精度的問題吧...看這題目,看這代碼,就沒有改的欲望...

hdu OJ 2449
http://acm.hdu.edu.cn/showproblem.php?pid=2449
哈爾濱現場賽的一道純高斯題,當時鶴牛敲了1個多小時...主要就是寫一個分數類,套個高精模板(偷懶點就Java...)搞定~~
注意下0和負數時的輸出即可。

fze OJ 1704
http://acm.fzu.edu.cn/problem.php?pid=1704
福大月賽的一道題目,還是經典的開關問題,但是方程數和變元數不同(考驗模板的時候到了~~),最后要求增廣陣的階,要用到高精度~~

Sgu 275 To xor or not to xor
http://acm.sgu.ru/problem.php?contest=0&problem=275
題解:
http://hi.baidu.com/czyuan%5Facm/blog/item/be3403d32549633d970a16ee.html

 

pps:對矩陣內涵的思考

 

如果不熟悉線性代數的概念,要去學習自然科學,現在看來就和文盲差不多。”,然而“按照現行的國際標准,線性代數是通過公理化來表述的,它是第二代數學模型,這就帶來了教學上的困難。”
* 矩陣究竟是什么東西?向量可以被認為是具有n個相互獨立的性質(維度)的對象的表示,矩陣又是什么呢?我們如果認為矩陣是一組列(行)向量組成的新的復合向量的展開式,那么為什么這種展開式具有如此廣泛的應用?特別是,為什么偏偏二維的展開式如此有用?如果矩陣中每一個元素又是一個向量,那么我們再展開一次,變成三維的立方陣,是不是更有用?
* 矩陣的乘法規則究竟為什么這樣規定?為什么這樣一種怪異的乘法規則卻能夠在實踐中發揮如此巨大的功效?很多看上去似乎是完全不相關的問題,最后竟然都歸結到矩陣的乘法,這難道不是很奇妙的事情?難道在矩陣乘法那看上去莫名其妙的規則下面,包含着世界的某些本質規律?如果是的話,這些本質規律是什么?
* 行列式究竟是一個什么東西?為什么會有如此怪異的計算規則?行列式與其對應方陣本質上是什么關系?為什么只有方陣才有對應的行列式,而一般矩陣就沒有(不要覺得這個問題很蠢,如果必要,針對m x n矩陣定義行列式不是做不到的,之所以不做,是因為沒有這個必要,但是為什么沒有這個必要)?而且,行列式的計算規則,看上去跟矩陣的任何計算規則都沒有直觀的聯系,為什么又在很多方面決定了矩陣的性質?難道這一切僅是巧合?
* 矩陣為什么可以分塊計算?分塊計算這件事情看上去是那么隨意,為什么竟是可行的?
* 對於矩陣轉置運算AT,有(AB)T = BTAT,對於矩陣求逆運算A-1,有(AB)-1 = B-1A-1。兩個看上去完全沒有什么關系的運算,為什么有着類似的性質?這僅僅是巧合嗎?
* 為什么說P-1AP得到的矩陣與A矩陣“相似”?這里的“相似”是什么意思?
* 特征值和特征向量的本質是什么?它們定義就讓人很驚訝,因為Ax =λx,一個諾大的矩陣的效應,竟然不過相當於一個小小的數λ,確實有點奇妙。但何至於用“特征”甚至“本征”來界定?它們刻划的究竟是什么?
今天先談談對線形空間和矩陣的幾個核心概念的理解。首先說說空間(space),這個概念是現代數學的命根子之一,從拓撲空間開始,一步步往上加定義,可以形成很多空間。線形空間其實還是比較初級的,如果在里面定義了范數,就成了賦范線性空間。賦范線性空間滿足完備性,就成了巴那赫空間;賦范線性空間中定義角度,就有了內積空間,內積空間再滿足完備性,就得到希爾伯特空間。
總之,空間有很多種。你要是去看某種空間的數學定義,大致都是“存在一個集合,在這個集合上定義某某概念,然后滿足某些性質”,就可以被稱為空間。這未免有點奇怪,為什么要用“空間”來稱呼一些這樣的集合呢?大家將會看到,其實這是很有道理的。
我們一般人最熟悉的空間,毫無疑問就是我們生活在其中的(按照牛頓的絕對時空觀)的三維空間,從數學上說,這是一個三維的歐幾里德空間,我們先不管那么多,先看看我們熟悉的這樣一個空間有些什么最基本的特點。仔細想想我們就會知道,這個三維的空間:1. 由很多(實際上是無窮多個)位置點組成;2. 這些點之間存在相對的關系;3. 可以在空間中定義長度、角度;4. 這個空間可以容納運動,這里我們所說的運動是從一個點到另一個點的移動(變換),而不是微積分意義上的“連續”性的運動,
事實上,不管是什么空間,都必須容納和支持在其中發生的符合規則的運動(變換)。你會發現,在某種空間中往往會存在一種相對應的變換,比如拓撲空間中有拓撲變換,線性空間中有線性變換,仿射空間中有仿射變換,其實這些變換都只不過是對應空間中允許的運動形式而已。
因此只要知道,“空間”是容納運動的一個對象集合,而變換則規定了對應空間的運動。
下面我們來看看線性空間。線性空間的定義任何一本書上都有,但是既然我們承認線性空間是個空間,那么有兩個最基本的問題必須首先得到解決,那就是:
1. 空間是一個對象集合,線性空間也是空間,所以也是一個對象集合。那么線性空間是什么樣的對象的集合?或者說,線性空間中的對象有什么共同點嗎?
2. 線性空間中的運動如何表述的?也就是,線性變換是如何表示的?
我們先來回答第一個問題,回答這個問題的時候其實是不用拐彎抹角的,可以直截了當的給出答案。線性空間中的任何一個對象,通過選取基和坐標的辦法,都可以表達為向量的形式。通常的向量空間我就不說了,舉兩個不那么平凡的例子:
L1. 最高次項不大於n次的多項式的全體構成一個線性空間,也就是說,這個線性空間中的每一個對象是一個多項式。如果我們以x0, x1, ..., xn為基,那么任何一個這樣的多項式都可以表達為一組n+1維向量,其中的每一個分量ai其實就是多項式中x(i-1)項的系數。值得說明的是,基的選取有多種辦法,只要所選取的那一組基線性無關就可以。這要用到后面提到的概念了,所以這里先不說,提一下而已。
L2. 閉區間[a, b]上的n階連續可微函數的全體,構成一個線性空間。也就是說,這個線性空間的每一個對象是一個連續函數。對於其中任何一個連續函數,根據魏爾斯特拉斯定理,一定可以找到最高次項不大於n的多項式函數,使之與該連續函數的差為0,也就是說,完全相等。這樣就把問題歸結為L1了。后面就不用再重復了。
所以說,向量是很厲害的,只要你找到合適的基,用向量可以表示線性空間里任何一個對象。這里頭大有文章,因為向量表面上只是一列數,但是其實由於它的有序性,所以除了這些數本身攜帶的信息之外,還可以在每個數的對應位置上攜帶信息。為什么在程序設計中數組最簡單,卻又威力無窮呢?根本原因就在於此。這是另一個問題了,這里就不說了。
下面來回答第二個問題,這個問題的回答會涉及到線性代數的一個最根本的問題。
線性空間中的運動,被稱為線性變換。也就是說,你從線性空間中的一個點運動到任意的另外一個點,都可以通過一個線性變化來完成。那么,線性變換如何表示呢?很有意思,在線性空間中,當你選定一組基之后,不僅可以用一個向量來描述空間中的任何一個對象,而且可以用矩陣來描述該空間中的任何一個運動(變換)。而使某個對象發生對應運動的方法,就是用代表那個運動的矩陣,乘以代表那個對象的向量。
簡而言之,在線性空間中選定基之后,向量刻畫對象,矩陣刻畫對象的運動,用矩陣與向量的乘法施加運動。
是的,矩陣的本質是運動的描述。如果以后有人問你矩陣是什么,那么你就可以響亮地告訴他,矩陣的本質是運動的描述
可是多么有意思啊,向量本身不是也可以看成是n x 1矩陣嗎?這實在是很奇妙,一個空間中的對象和運動竟然可以用相類同的方式表示。能說這是巧合嗎?如果是巧合的話,那可真是幸運的巧合!可以說,線性代數中大多數奇妙的性質,均與這個巧合有直接的關系。
   “矩陣是運動的描述”,到現在為止,好像大家都還沒什么意見。但是我相信早晚會有數學系出身的網友來拍板轉。因為運動這個概念,在數學和物理里是跟微積分聯系在一起的。我們學習微積分的時候,總會有人照本宣科地告訴你,初等數學是研究常量的數學,是研究靜態的數學,高等數學是變量的數學,是研究運動的數學。大家口口相傳,差不多人人都知道這句話。但是真知道這句話說的是什么意思的人,好像也不多。簡而言之,在我們人類的經驗里,運動是一個連續過程,從A點到B點,就算走得最快的光,也是需要一個時間來逐點地經過AB之間的路徑,這就帶來了連續性的概念。而連續這個事情,如果不定義極限的概念,根本就解釋不了。古希臘人的數學非常強,但就是缺乏極限觀念,所以解釋不了運動,被芝諾的那些著名悖論(飛箭不動、飛毛腿阿喀琉斯跑不過烏龜等四個悖論)搞得死去活來。因為這篇文章不是講微積分的,所以我就不多說了。有興趣的讀者可以去看看齊民友教授寫的《重溫微積分》。我就是讀了這本書開頭的部分,才明白“高等數學是研究運動的數學”這句話的道理。
不過在我這個《理解矩陣》的文章里,“運動”的概念不是微積分中的連續性的運動,而是瞬間發生的變化。比如這個時刻在A點,經過一個“運動”,一下子就“躍遷”到了B點,其中不需要經過A點與B點之間的任何一個點。這樣的“運動”,或者說“躍遷”,是違反我們日常的經驗的。不過了解一點量子物理常識的人,就會立刻指出,量子(例如電子)在不同的能量級軌道上跳躍,就是瞬間發生的,具有這樣一種躍遷行為。所以說,自然界中並不是沒有這種運動現象,只不過宏觀上我們觀察不到。但是不管怎么說,“運動”這個詞用在這里,還是容易產生歧義的,說得更確切些,應該是“躍遷”。因此這句話可以改成:
“矩陣是線性空間里躍遷的描述”。
可是這樣說又太物理,也就是說太具體,而不夠數學,也就是說不夠抽象。因此我們最后換用一個正牌的數學術語——變換,來描述這個事情。這樣一說,大家就應該明白了,所謂變換,其實就是空間里從一個點(元素/對象)到另一個點(元素/對象)的躍遷。比如說,拓撲變換,就是在拓撲空間里從一個點到另一個點的躍遷。再比如說,仿射變換,就是在仿射空間里從一個點到另一個點的躍遷。附帶說一下,這個仿射空間跟向量空間是親兄弟。做計算機圖形學的朋友都知道,盡管描述一個三維對象只需要三維向量,但所有的計算機圖形學變換矩陣都是4 x 4的。說其原因,很多書上都寫着“為了使用中方便”,這在我看來簡直就是企圖蒙混過關。真正的原因,是因為在計算機圖形學里應用的圖形變換,實際上是在仿射空間而不是向量空間中進行的。想想看,在向量空間里相一個向量平行移動以后仍是相同的那個向量,而現實世界等長的兩個平行線段當然不能被認為同一個東西,所以計算機圖形學的生存空間實際上是仿射空間。而仿射變換的矩陣表示根本就是4 x 4的。又扯遠了,有興趣的讀者可以去看《計算機圖形學——幾何工具算法詳解》。
一旦我們理解了“變換”這個概念,矩陣的定義就變成:
“矩陣是線性空間里的變換的描述。”
到這里為止,我們終於得到了一個看上去比較數學的定義。不過還要多說幾句。教材上一般是這么說的,在一個線性空間V里的一個線性變換T,當選定一組基之后,就可以表示為矩陣。因此我們還要說清楚到底什么是線性變換,什么是基,什么叫選定一組基。線性變換的定義是很簡單的,設有一種變換T,使得對於線性空間V中間任何兩個不相同的對象x和y,以及任意實數a和b,有:
T(ax + by) = aT(x) + bT(y),
那么就稱T為線性變換。
定義都是這么寫的,但是光看定義還得不到直覺的理解。線性變換究竟是一種什么樣的變換?我們剛才說了,變換是從空間的一個點躍遷到另一個點,而線性變換,就是從一個線性空間V的某一個點躍遷到另一個線性空間W的另一個點的運動。這句話里蘊含着一層意思,就是說一個點不僅可以變換到同一個線性空間中的另一個點,而且可以變換到另一個線性空間中的另一個點去。不管你怎么變,只要變換前后都是線性空間中的對象,這個變換就一定是線性變換,也就一定可以用一個非奇異矩陣來描述。而你用一個非奇異矩陣去描述的一個變換,一定是一個線性變換。有的人可能要問,這里為什么要強調非奇異矩陣?所謂非奇異,只對方陣有意義,那么非方陣的情況怎么樣?這個說起來就會比較冗長了,最后要把線性變換作為一種映射,並且討論其映射性質,以及線性變換的核與像等概念才能徹底講清楚。我覺得這個不算是重點,如果確實有時間的話,以后寫一點。以下我們只探討最常用、最有用的一種變換,就是在同一個線性空間之內的線性變換。也就是說,下面所說的矩陣,不作說明的話,就是方陣,而且是非奇異方陣。學習一門學問,最重要的是把握主干內容,迅速建立對於這門學問的整體概念,不必一開始就考慮所有的細枝末節和特殊情況,自亂陣腳。
接着往下說,什么是基呢?這個問題在后面還要大講一番,這里只要把基看成是線性空間里的坐標系就可以了。注意是坐標系,不是坐標值,這兩者可是一個“對立矛盾統一體”。這樣一來,“選定一組基”就是說在線性空間里選定一個坐標系。就這意思。
好,最后我們把矩陣的定義完善如下:
“矩陣是線性空間中的線性變換的一個描述。在一個線性空間中,只要我們選定一組基,那么對於任何一個線性變換,都能夠用一個確定的矩陣來加以描述。”
理解這句話的關鍵,在於把“線性變換”與“線性變換的一個描述”區別開。一個是那個對象,一個是對那個對象的表述。就好像我們熟悉的面向對象編程中,一個對象可以有多個引用,每個引用可以叫不同的名字,但都是指的同一個對象。如果還不形象,那就干脆來個很俗的類比。
比如有一頭豬,你打算給它拍照片,只要你給照相機選定了一個鏡頭位置,那么就可以給這頭豬拍一張照片。這個照片可以看成是這頭豬的一個描述,但只是一個片面的的描述,因為換一個鏡頭位置給這頭豬拍照,能得到一張不同的照片,也是這頭豬的另一個片面的描述。所有這樣照出來的照片都是這同一頭豬的描述,但是又都不是這頭豬本身。
同樣的,對於一個線性變換,只要你選定一組基,那么就可以找到一個矩陣來描述這個線性變換。換一組基,就得到一個不同的矩陣。所有這些矩陣都是這同一個線性變換的描述,但又都不是線性變換本身。
但是這樣的話,問題就來了如果你給我兩張豬的照片,我怎么知道這兩張照片上的是同一頭豬呢?同樣的,你給我兩個矩陣,我怎么知道這兩個矩陣是描述的同一個線性變換呢?如果是同一個線性變換的不同的矩陣描述,那就是本家兄弟了,見面不認識,豈不成了笑話。
好在,我們可以找到同一個線性變換的矩陣兄弟們的一個性質,那就是:
若矩陣A與B是同一個線性變換的兩個不同的描述(之所以會不同,是因為選定了不同的基,也就是選定了不同的坐標系),則一定能找到一個非奇異矩陣P,使得A、B之間滿足這樣的關系:
A = P-1BP
線性代數稍微熟一點的讀者一下就看出來,這就是相似矩陣的定義。沒錯,所謂相似矩陣,就是同一個線性變換的不同的描述矩陣。按照這個定義,同一頭豬的不同角度的照片也可以成為相似照片。俗了一點,不過能讓人明白。
而在上面式子里那個矩陣P,其實就是A矩陣所基於的基與B矩陣所基於的基這兩組基之間的一個變換關系。關於這個結論,可以用一種非常直覺的方法來證明(而不是一般教科書上那種形式上的證明),如果有時間的話,我以后在blog里補充這個證明。
這個發現太重要了。原來一族相似矩陣都是同一個線性變換的描述啊!難怪這么重要!工科研究生課程中有矩陣論、矩陣分析等課程,其中講了各種各樣的相似變換,比如什么相似標准型,對角化之類的內容,都要求變換以后得到的那個矩陣與先前的那個矩陣式相似的,為什么這么要求?因為只有這樣要求,才能保證變換前后的兩個矩陣是描述同一個線性變換的。當然,同一個線性變換的不同矩陣描述,從實際運算性質來看並不是不分好環的。有些描述矩陣就比其他的矩陣性質好得多。這很容易理解,同一頭豬的照片也有美丑之分嘛。所以矩陣的相似變換可以把一個比較丑的矩陣變成一個比較美的矩陣,而保證這兩個矩陣都是描述了同一個線性變換。
這樣一來,矩陣作為線性變換描述的一面,基本上說清楚了。但是,事情沒有那么簡單,或者說,線性代數還有比這更奇妙的性質,那就是,矩陣不僅可以作為線性變換的描述,而且可以作為一組基的描述。而作為變換的矩陣,不但可以把線性空間中的一個點給變換到另一個點去,而且也能夠把線性空間中的一個坐標系(基)表換到另一個坐標系(基)去。而且,變換點與變換坐標系,具有異曲同工的效果。線性代數里最有趣的奧妙,就蘊含在其中。理解了這些內容,線性代數里很多定理和規則會變得更加清晰、直覺。
首先來總結一下前面兩部分的一些主要結論:
1. 首先有空間,空間可以容納對象運動的。一種空間對應一類對象。
2. 有一種空間叫線性空間,線性空間是容納向量對象運動的。
3. 運動是瞬時的,因此也被稱為變換。
4. 矩陣是線性空間中運動(變換)的描述。
5. 矩陣與向量相乘,就是實施運動(變換)的過程。
6. 同一個變換,在不同的坐標系下表現為不同的矩陣,但是它們的本質是一樣的,所以本征值相同。
   下面讓我們把視力集中到一點以改變我們以往看待矩陣的方式。我們知道,線性空間里的基本對象是向量,而向量是這么表示的:
        [a1, a2, a3, ..., an]
 矩陣呢?矩陣是這么表示的:
        a11, a12, a13, ..., a1n
        a21, a22, a23, ..., a2n
                  ...
        an1, an2, an3, ..., ann
        不用太聰明,我們就能看出來,矩陣是一組向量組成的。特別的,n維線性空間里的方陣是由n個n維向量組成的。我們在這里只討論這個n階的、非奇異的方陣,如果一組向量是彼此線性無關的話,那么它們就可以成為度量這個線性空間的一組基,從而事實上成為一個坐標系體系,其中每一個向量都躺在一根坐標軸上,並且成為那根坐標軸上的基本度量單位(長度1)。現在到了關鍵的一步。看上去矩陣就是由一組向量組成的,而且如果矩陣非奇異的話(我說了,只考慮這種情況),那么組成這個矩陣的那一組向量也就是線性無關的了,也就可以成為度量線性空間的一個坐標系。結論:矩陣描述了一個坐標系。之所以矩陣又是運動,又是坐標系,那是因為——“運動等價於坐標系變換”。對不起,這話其實不准確,我只是想讓你印象深刻。准確的說法是:“對象的變換等價於坐標系的變換”。或者:“固定坐標系下一個對象的變換等價於固定對象所處的坐標系變換。” 說白了就是: “運動是相對的。”        
    讓我們想想,達成同一個變換的結果,比如把點(1, 1)變到點(2, 3)去,你可以有兩種做法。第一,坐標系不動,點動,把(1, 1)點挪到(2, 3)去。第二,點不動,變坐標系,讓x軸的度量(單位向量)變成原來的1/2,讓y軸的度量(單位向量)變成原先的1/3,這樣點還是那個點,可是點的坐標就變成(2, 3)了。方式不同,結果一樣。從第一個方式來看,那就是我在《理解矩陣》1/2中說的,把矩陣看成是運動描述,矩陣與向量相乘就是使向量(點)運動的過程。在這個方式下,
      Ma = b的意思是:
      “向量a經過矩陣M所描述的變換,變成了向量b。”而從第二個方式來看,矩陣M描述了一個坐標系,姑且也稱之為M。那么:
      Ma = b的意思是:
      “有一個向量,它在坐標系M的度量下得到的度量結果向量為a,那么它在坐標系I的度量下,這個向量的度量結果是b。”
      這里的I是指單位矩陣,就是主對角線是1,其他為零的矩陣。而這兩個方式本質上是等價的。我希望你務必理解這一點,因為這是本篇的關鍵。正因為是關鍵,所以我得再解釋一下。在M為坐標系的意義下,如果把M放在一個向量a的前面,形成Ma的樣式,我們可以認為這是對向量a的一個環境聲明。它相當於是說: “注意了!這里有一個向量,它在坐標系M中度量,得到的度量結果可以表達為a。可是它在別的坐標系里度量的話,就會得到不同的結果。為了明確,我把M放在前面,讓你明白,這是該向量在坐標系M中度量的結果。” 那么我們再看孤零零的向量b:
       b          多看幾遍,你沒看出來嗎?它其實不是b,它是:
       Ib
     也就是說:“在單位坐標系,也就是我們通常說的直角坐標系I中,有一個向量,度量的結果是b。”
       而 Ma = Ib的意思就是說:
     “在M坐標系里量出來的向量a,跟在I坐標系里量出來的向量b,其實根本就是一個向量啊!”這哪里是什么乘法計算,根本就是身份識別嘛。從這個意義上我們重新理解一下向量。向量這個東西客觀存在,但是要把它表示出來,就要把它放在一個坐標系中去度量它,然后把度量的結果(向量在各個坐標軸上的投影值)按一定順序列在一起,就成了我們平時所見的向量表示形式。你選擇的坐標系(基)不同,得出來的向量的表示就不同。向量還是那個向量,選擇的坐標系不同,其表示方式就不同。因此,按道理來說,每寫出一個向量的表示,都應該聲明一下這個表示是在哪個坐標系中度量出來的。表示的方式,就是 Ma,也就是說,有一個向量,在M矩陣表示的坐標系中度量出來的結果為a。我們平時說一個向量是[2 3 5 7]T,隱含着是說,這個向量在 I 坐標系中的度量結果是[2 3 5 7]T,因此,這個形式反而是一種簡化了的特殊情況。

注意到,M矩陣表示出來的那個坐標系,由一組基組成,而那組基也是由向量組成的,同樣存在這組向量是在哪個坐標系下度量而成的問題。也就是說,表述一個矩陣的一般方法,也應該要指明其所處的基准坐標系。所謂M,其實是 IM,也就是說,M中那組基的度量是在 I 坐標系中得出的。從這個視角來看,M×N也不是什么矩陣乘法了,而是聲明了一個在M坐標系中量出的另一個坐標系N,其中M本身是在I坐標系中度量出來的。
回過頭來說變換的問題。我剛才說,“固定坐標系下一個對象的變換等價於固定對象所處的坐標系變換”,那個“固定對象”我們找到了,就是那個向量。但是坐標系的變換呢?我怎么沒看見?
請看:
Ma = Ib
我現在要變M為I,怎么變?對了,再前面乘以個M-1,也就是M的逆矩陣。換句話說,你不是有一個坐標系M嗎,現在我讓它乘以個M-1,變成I,這樣一來的話,原來M坐標系中的a在I中一量,就得到b了。
我建議你此時此刻拿起紙筆,畫畫圖,求得對這件事情的理解。比如,你畫一個坐標系,x軸上的衡量單位是2,y軸上的衡量單位是3,在這樣一個坐標系里,坐標為(1,1)的那一點,實際上就是笛卡爾坐標系里的點(2, 3)。而讓它原形畢露的辦法,就是把原來那個坐標系:
2 0
0 3
的x方向度量縮小為原來的1/2,而y方向度量縮小為原來的1/3,這樣一來坐標系就變成單位坐標系I了。保持點不變,那個向量現在就變成了(2, 3)了。
怎么能夠讓“x方向度量縮小為原來的1/2,而y方向度量縮小為原來的1/3”呢?就是讓原坐標系:
2 0
0 3
被矩陣:
1/2 0
0 1/3
左乘。而這個矩陣就是原矩陣的逆矩陣。
下面我們得出一個重要的結論:
“對坐標系施加變換的方法,就是讓表示那個坐標系的矩陣與表示那個變化的矩陣相乘。”
再一次的,矩陣的乘法變成了運動的施加。只不過,被施加運動的不再是向量,而是另一個坐標系。
如果你覺得你還搞得清楚,請再想一下剛才已經提到的結論,矩陣MxN,一方面表明坐標系N在運動M下的變換結果,另一方面,把M當成N的前綴,當成N的環境描述,那么就是說,在M坐標系度量下,有另一個坐標系N。這個坐標系N如果放在I坐標系中度量,其結果為坐標系MxN。
在這里,我實際上已經回答了一般人在學習線性代數是最困惑的一個問題,那就是為什么矩陣的乘法要規定成這樣。簡單地說,是因為:
1. 從變換的觀點看,對坐標系N施加M變換,就是把組成坐標系N的每一個向量施加M變換。
2. 從坐標系的觀點看,在M坐標系中表現為N的另一個坐標系,這也歸結為,對N坐標系基的每一個向量,把它在I坐標系中的坐標找出來,然后匯成一個新的矩陣。
3. 至於矩陣乘以向量為什么要那樣規定,那是因為一個在M中度量為a的向量,如果想要恢復在I中的真像,就必須分別與M中的每一個向量進行內積運算。我把這個結論的推導留給感興趣的朋友吧。應該說,其實到了這一步,已經很容易了。
綜合以上1/2/3,矩陣的乘法就得那么規定,一切有根有據,絕不是哪個神經病胡思亂想出來的。

我已經無法說得更多了。矩陣又是坐標系,又是變換。到底是坐標系,還是變換,已經說不清楚了,運動與實體在這里統一了,物質與意識的界限已經消失了,一切歸於無法言說,無法定義了。道可道,非常道,名可名,非常名。矩陣是在是不可道之道,不可名之名的東西。到了這個時候,我們不得不承認,我們偉大的線性代數課本上說的矩陣定義,是無比正確的:
“矩陣就是由m行n列數放在一起組成的數學對象。”
好了,這基本上就是我想說的全部了。還留下一個行列式的問題。矩陣M的行列式實際上是組成M的各個向量按照平行四邊形法則搭成一個n維立方體的體積。對於這一點,我只能感嘆於其精妙,卻無法揭開其中奧秘了。也許我掌握的數學工具不夠,我希望有人能夠給我們大家講解其中的道理了。
我不知道是否講得足夠清楚了,反正這一部分需要您花些功夫去推敲。
此外,請大家不必等待這個系列的后續部分。以我的工作情況而言,近期內很難保證繼續投入腦力到這個領域中,盡管我仍然對此興致濃厚。不過如果還有(四)的話,可能是一些站在應用層面的考慮,比如對計算機圖形學相關算法的理解。但是我不承諾這些討論近期內會出現了。

 

“分”的反義字是“和”,是我們熟悉的字。比如:2+3=5,從左往右運算,我們叫求和。那么“分”呢,既然是反義字,就把上面的等式反過來:5=2+3。
把一個對象表示成兩個以至更多的對象的和,這個過程叫分析。
通常來說,分析對象應當與被分析對象一致。是數就都是數,是函數就都是函數,是向量就都是向量,是矩陣就都是矩陣。

求和是數學里最基本的運算,減、乘、除是從求和中衍生出來的。而更高級的冪、指、對、三角、微積分等,也是一層一層建立起來的, 最根本的還是這個求和。求和最簡單,最容易計算,性質也最簡單。所以成了分析的基本出發點。

分析的妙處在於,通過分析可以將較復雜的對象划分為較簡單的對象。比如2和3就比5簡單。單獨研究2的性質,再單獨研究3的性質,再通過簡單的求和,就可以把握5的性質。把復雜的東西划分成若干簡單對象的和,對各簡單對象搞各個擊破,再加起來,復雜的東西也就被掌握了。

分析是西方思想中一個根本性的東西。西方人認為,事物總是有因果的,看到了結果,要分析原因。所謂分析原因,就是找出一堆因素,說明這堆因素合起來導致了結果。西方人認為,事物總是可以分析的。看到了整體,就要把那些合成這個整體的局部一一分析出來。 現代科學很大一部分就是這么回事。

大學數學里,有很多內容就是在講分析。數學里的分析還要把含義拓展,就是把一個數學對象合理地表示成若干更簡單對象與實數系數之積的和。但微積分和線性代數各有側重。微積分研究的是無窮項求和。無窮項之和與有窮項之和是本質不同的。但是無窮項之和是無法運算的,至少不實際。所以要想辦法通過一種辦法用有窮項之和來近似的代替,這就是逼近。逼近成立的條件是收斂,就是說,只有從一個收斂的無窮項的開頭截出一部分來求和,才能被認為是逼近。華人數學家項武義說,微積分就逼近這一板斧,但是無往而不利。
微積分主要研究函數,連續函數的因變量y會由於自變量x的變化而變化。這種變化也是要分析的。當x從 x0變成x1時,y是怎樣從y0變到y1 的?按照上面的說法,“y的變化(y1-y0)”這一個數學對象,要用一系列比較簡單的“變化”相加來表示。數學家找到了一個收斂的“變化”對象的序列,排在頭一位的是一個線性的變化量,它的系數就是導數,它本身就是微分dy。數學家又發現,當x的變化量無窮小時,從這個無窮的、收斂的“變化”對象序列中,只要截出第一項,也就是微分dy,就無論如何可以精確描述y的變化了。曾在一本書上見過這樣的說法,泰勒公式是數學分析的頂峰。不知道是不是有道理。我自己覺得是這么回事。有了泰勒公式,我們可以任意精確地算一個函數在某一點上的值。畢竟只是實數求和嘛。

但是為了表示泰勒公式,我們卻用了一個挺復雜的連加代數式。代數式不能象實數那樣簡單加起來得到一個對象,它只能表示成和的形式。這是我們意識到,在這個連加式中各對象存在某些特別的不同,使它們沒法簡單地加到一起。 因此我們有必要討論,把一些性質不同的東西加到一起所形成的這個對象有什么性質。 這就是向量。

微積分研究如何把一個對象分解為無窮項同質對象之和,線性代數研究“有限項異質對象之和”這個新對象的性質。一方面,上面說過,微積分到最后還是要化無窮為有窮,化精確為逼近;另一方面,異質對象經過某種處理可以轉化為同質對象。比如不同次的冪函數是異質對象,但是一旦代入具體數值則都可以轉化為實數,變成了同質對象。因此線性代數研究的問題對微積分很重要。故我認為大學里應先講線性代數,后講微積分。

我們的微積分教學,將重點過分傾注在微分和積分的運算上了,其實實踐中更為重要的是我們稱為“級數”的那部分內容。即研究如何將一個量表達為一個數項級數,如何將一個函數表達為一個函數項級數。

線性代數把異質對象之和(向量)作為研究的基礎,研究這些新定義的對象加起來又可以表示什么。其結論是,有限數量的向量連加起來,有可能具有這樣的能力,即同維的全部向量都可以表示成這些向量的和。這樣的一組具有充分表現能力的向量,是線性無關的向量,組成了一個向量空間,而它們自己構成了這個向量空間里的一組基。

回到分析的概念上,一個向量總可以表示為若干個同階向量之和,這就是向量的分析。但是並不是所有的這些分析都具有相同的價值。在某種運算中,某種特別的分析能夠提供特別優越的性,從而大大簡化運算。比如在大多數情況下,將一個向量表示成一組單位正交基向量的和,就能夠在計算中獲得特別的便利。面對某個問題,尋找一個最優越的分析形式,把要研究的對象合理地表示成具有特殊性質的基對象與實數系數之積的和,這是分析的重要步驟,也是成功的關鍵。在這種表示式中,系數稱為坐標。

經典的方法都是以找到一組性質優良的基為開端的,例如:
傅立葉分析以正交函數系為基,因此具有優良性質,自1904年以來取代冪函數系,成為分析主流。
在曲線和曲面擬合中,正交多項式集構成了最佳基函數。 拉格朗日插值多項式具有一個特別的性質,即在本結點上為1,在其他結點上為0。
有限元中的形函數類似拉氏插值多項式。
結構動力學中的主振型迭加法,也是以相互正交的主振型為基,對多質點體系位移進行分析的。

舉兩個例子:說到采樣,大家的第一反應肯定是一個詞“2倍”(采樣定理)。學得比較扎實的,可能還會把為什么是2倍解釋清楚。但我對采樣的理解是:采樣實際上是在進行正交分解,采樣值不過是在一組正交基下分解的系數。如果原信號屬於該組正交基所張成的線性子空間,那么該信號就能無失真的恢復(滿足采樣定理)。學過信號處理的朋友,你知道這組正交基是什么嗎?:)第二個例子是關於為什么傅里葉變換在線性系統理論中如此重要?答案可能五花八門,但我認為我的理解是比較深入的:原因是傅里葉基是所有線性時不變算子的特征向量(和本文聯系起來了)。這句話解釋起來比較費工夫,但是傅里葉變換能和特征向量聯系起來,大家一定感覺很有趣吧。

 

 

特征向量確實有很明確的幾何意義,矩陣(既然討論特征向量的問題,當然是方陣,這里不討論廣義特征向量的概念,就是一般的特征向量)乘以一個向量的結果仍是同維數的一個向量,因此,矩陣乘法對應了一個變換,把一個向量變成同維數的另一個向量,那么變換的效果是什么呢?這當然與方陣的構造有密切關系,比如可以取適當的二維方陣,使得這個變換的效果就是將平面上的二維向量逆時針旋轉30度,這時我們可以問一個問題,有沒有向量在這個變換下不改變方向呢?可以想一下,除了零向量,沒有其他向量可以在平面上旋轉30度而不改變方向的,所以這個變換對應的矩陣(或者說這個變換自身)沒有特征向量(注意:特征向量不能是零向量),所以一個變換的特征向量是這樣一種向量,它經過這種特定的變換后保持方向不變,只是進行長度上的伸縮而已(再想想特征向量的原始定義Ax=cx,你就恍然大悟了,看到了嗎?cx是方陣A對向量x進行變換后的結果,但顯然cx和x的方向相同),而且x是特征向量的話,ax也是特征向量(a是標量且不為零),所以所謂的特征向量不是一個向量而是一個向量族,另外,特征值只不過反映了特征向量在變換時的伸縮倍數而已,對一個變換而言,特征向量指明的方向才是很重要的,特征值不是那么重要,雖然我們求這兩個量時先求出特征值,但特征向量才是更本質的東西!

比如平面上的一個變換,把一個向量關於橫軸做鏡像對稱變換,即保持一個向量的橫坐標不變,但縱坐標取相反數,把這個變換表示為矩陣就是[1 0;0 -1],其中分號表示換行,顯然[1 0;0 -1]*[a b]'=[a -b]',其中上標'表示取轉置,這正是我們想要的效果,那么現在可以猜一下了,這個矩陣的特征向量是什么?想想什么向量在這個變換下保持方向不變,顯然,橫軸上的向量在這個變換下保持方向不變(記住這個變換是鏡像對稱變換,那鏡子表面上(橫軸上)的向量當然不會變化),所以可以直接猜測其特征向量是 [a 0]'(a不為0),還有其他的嗎?有,那就是縱軸上的向量,這時經過變換后,其方向反向,但仍在同一條軸上,所以也被認為是方向沒有變化,所以[0 b]'(b不為0)也是其特征向量,去求求矩陣[1 0;0 -1]的特征向量就知道對不對了!

//最近一直在做一個數論專題,后期有待整理,先將有關資料收藏下,在學習高斯消元的時候看了czyuan大牛的此文獲益匪淺,czyuan的此份模板可以解決大多高斯問題,當然並不是萬能的,其中建立矩陣是難點,需要自己琢磨,並且對於方程組是否有解、解的個數以及自由元等問題也需要自己做題慢慢思考,自己做了兩三道題前前后后在建矩陣以及對一些解的問題在Gauss函數中改了幾十次,逐漸摸索,還不算掌握的好,有待再慢慢練。

 

高斯消元法,是線性代數中的一個算法,可用來求解線性方程組,並可以求出矩陣的秩,以及求出可逆方陣的逆矩陣。
高斯消元法的原理是:
若用初等行變換將增廣矩陣 化為 ,則AX = B與CX = D是同解方程組。
所以我們可以用初等行變換把增廣矩陣轉換為行階梯陣,然后回代求出方程的解。

以上是線性代數課的回顧,下面來說說高斯消元法在編程中的應用。

首先,先介紹程序中高斯消元法的步驟:
(我們設方程組中方程的個數為equ,變元的個數為var,注意:一般情況下是n個方程,n個變元,但是有些題目就故意讓方程數與變元數不同)

1. 把方程組轉換成增廣矩陣。

2. 利用初等行變換來把增廣矩陣轉換成行階梯陣。
枚舉k從0到equ – 1,當前處理的列為col(初始為0) ,每次找第k行以下(包括第k行),col列中元素絕對值最大的列與第k行交換。如果col列中的元素全為0,那么則處理col + 1列,k不變。

3. 轉換為行階梯陣,判斷解的情況。

① 無解
當方程中出現(0, 0, …, 0, a)的形式,且a != 0時,說明是無解的。

② 唯一解
條件是k = equ,即行階梯陣形成了嚴格的上三角陣。利用回代逐一求出解集。

③ 無窮解。
條件是k < equ,即不能形成嚴格的上三角形,自由變元的個數即為equ – k,但有些題目要求判斷哪些變元是不缺定的。
    這里單獨介紹下這種解法:
首先,自由變元有var - k個,即不確定的變元至少有var - k個。我們先把所有的變元視為不確定的。在每個方程中判斷不確定變元的個數,如果大於1個,則該方程無法求解。如果只有1個變元,那么該變元即可求出,即為確定變元。

以上介紹的是求解整數線性方程組的求法,復雜度是O(n3)。浮點數線性方程組的求法類似,但是要在判斷是否為0時,加入EPS,以消除精度問題。

下面講解幾道OJ上的高斯消元法求解線性方程組的題目:

POJ 1222 EXTENDED LIGHTS OUT
http://acm.pku.edu.cn/JudgeOnline/problem?id=1222
POJ 1681 Painter's Problem
http://acm.pku.edu.cn/JudgeOnline/problem?id=1681
POJ 1753 Flip Game
http://acm.pku.edu.cn/JudgeOnline/problem?id=1753
POJ 1830 開關問題
http://acm.pku.edu.cn/JudgeOnline/problem?id=1830

POJ 3185 The Water Bowls

http://acm.pku.edu.cn/JudgeOnline/problem?id=3185
開關窗戶,開關燈問題,很典型的求解線性方程組的問題。方程數和變量數均為行數*列數,直接套模板求解即可。但是,當出現無窮解時,需要枚舉解的情況,因為無法判斷哪種解是題目要求最優的。

POJ 2947 Widget Factory
http://acm.pku.edu.cn/JudgeOnline/problem?id=2947
求解同余方程組問題。與一般求解線性方程組的問題類似,只要在求解過程中加入取余即可。
注意:當方程組唯一解時,求解過程中要保證解在[3, 9]之間。

POJ 1166 The Clocks
http://acm.pku.edu.cn/JudgeOnline/problem?id=1166
經典的BFS問題,有各種解法,也可以用逆矩陣進行矩陣相乘。
但是這道題用高斯消元法解決好像有些問題(困擾了我N天...持續困擾中...),由於周期4不是素數,故在求解過程中不能進行取余(因為取余可能導致解集變大),但最后求解集時,還是需要進行取余操作,那么就不能保證最后求出的解是正確的...在discuss里提問了好幾天也沒人回答...希望哪位路過的大牛指點下~~

POJ 2065 SETI
http://acm.pku.edu.cn/JudgeOnline/problem?id=2065
同樣是求解同余方程組問題,由於題目中的p是素數,可以直接在求解時取余,套用模板求解即可。(雖然AC的人很少,但它還是比較水的一道題,)

POJ 1487 Single-Player Games
http://acm.pku.edu.cn/JudgeOnline/problem?id=1487
很麻煩的一道題目...題目中的敘述貌似用到了編譯原理中的詞法定義(看了就給人不想做的感覺...)
解方程組的思想還是很好看出來了(前提是通讀題目不下5遍...),但如果把樹的字符串表達式轉換成方程組是個難點,我是用棧 + 遞歸的做法分解的。首先用棧的思想求出該結點的孩子數,然后遞歸分別求解各個孩子。
這題解方程組也與眾不同...首先是求解浮點數方程組,要注意精度問題,然后又詢問不確定的變元,按前面說的方法求解。
一頓折騰后,這題居然寫了6000+B...而且囧的是巨人C++ WA,G++ AC,可能還是精度的問題吧...看這題目,看這代碼,就沒有改的欲望...

hdu OJ 2449
http://acm.hdu.edu.cn/showproblem.php?pid=2449
哈爾濱現場賽的一道純高斯題,當時鶴牛敲了1個多小時...主要就是寫一個分數類,套個高精模板(偷懶點就Java...)搞定~~
注意下0和負數時的輸出即可。

fze OJ 1704
http://acm.fzu.edu.cn/problem.php?pid=1704
福大月賽的一道題目,還是經典的開關問題,但是方程數和變元數不同(考驗模板的時候到了~~),最后要求增廣陣的階,要用到高精度~~

Sgu 275 To xor or not to xor
http://acm.sgu.ru/problem.php?contest=0&problem=275
題解:
http://hi.baidu.com/czyuan%5Facm/blog/item/be3403d32549633d970a16ee.html

這里提供下自己寫的還算滿意的求解整數線性方程組的模板(浮點數類似,就不提供了)~~

  1 /* 用於求整數解得方程組. */
  2 
  3 #include <iostream>
  4 #include <string>
  5 #include <cmath>
  6 using namespace std;
  7 
  8 const int maxn = 105;
  9 
 10 int equ, var; // 有equ個方程,var個變元。增廣陣行數為equ, 分別為0到equ - 1,列數為var + 1,分別為0到var.
 11 int a[maxn][maxn];
 12 int x[maxn]; // 解集.
 13 bool free_x[maxn]; // 判斷是否是不確定的變元.
 14 int free_num;
 15 
 16 void Debug(void)
 17 {
 18     int i, j;
 19     for (i = 0; i < equ; i++)
 20     {
 21         for (j = 0; j < var + 1; j++)
 22         {
 23             cout << a[i][j] << " ";
 24         }
 25         cout << endl;
 26     }
 27     cout << endl;
 28 }
 29 
 30 inline int gcd(int a, int b)
 31 {
 32     int t;
 33     while (b != 0)
 34     {
 35         t = b;
 36         b = a % b;
 37         a = t;
 38     }
 39     return a;
 40 }
 41 
 42 inline int lcm(int a, int b)
 43 {
 44     return a * b / gcd(a, b);
 45 }
 46 
 47 // 高斯消元法解方程組(Gauss-Jordan elimination).(-2表示有浮點數解,但無整數解,-1表示無解,0表示唯一解,大於0表示無窮解,並返回自由變元的個數)
 48 int Gauss(void)
 49 {
 50     int i, j, k;
 51     int max_r; // 當前這列絕對值最大的行.
 52 int col; // 當前處理的列.
 53     int ta, tb;
 54     int LCM;
 55     int temp;
 56     int free_x_num;
 57     int free_index;
 58     // 轉換為階梯陣.
 59     col = 0; // 當前處理的列.
 60     for (k = 0; k < equ && col < var; k++, col++)
 61     { // 枚舉當前處理的行.
 62         // 找到該col列元素絕對值最大的那行與第k行交換.(為了在除法時減小誤差)
 63         max_r = k;
 64         for (i = k + 1; i < equ; i++)
 65         {
 66             if (abs(a[i][col]) > abs(a[max_r][col])) max_r = i;
 67         }
 68         if (max_r != k)
 69         { // 與第k行交換.
 70             for (j = k; j < var + 1; j++) swap(a[k][j], a[max_r][j]);
 71         }
 72         if (a[k][col] == 0)
 73         { // 說明該col列第k行以下全是0了,則處理當前行的下一列.
 74             k--; continue;
 75         }
 76         for (i = k + 1; i < equ; i++)
 77         { // 枚舉要刪去的行.
 78             if (a[i][col] != 0)
 79     {
 80                 LCM = lcm(abs(a[i][col]), abs(a[k][col]));
 81                 ta = LCM / abs(a[i][col]), tb = LCM / abs(a[k][col]);
 82                 if (a[i][col] * a[k][col] < 0) tb = -tb; // 異號的情況是兩個數相加.
 83                 for (j = col; j < var + 1; j++)
 84                 {
 85                     a[i][j] = a[i][j] * ta - a[k][j] * tb;
 86                 }
 87     }
 88         }
 89     }
 90     Debug();
 91     // 1. 無解的情況: 化簡的增廣陣中存在(0, 0, ..., a)這樣的行(a != 0).
 92     for (i = k; i < equ; i++)
 93     { // 對於無窮解來說,如果要判斷哪些是自由變元,那么初等行變換中的交換就會影響,則要記錄交換.
 94         if (a[i][col] != 0) return -1;
 95     }
 96     // 2. 無窮解的情況: 在var * (var + 1)的增廣陣中出現(0, 0, ..., 0)這樣的行,即說明沒有形成嚴格的上三角陣.
 97     // 且出現的行數即為自由變元的個數.
 98     if (k < var)
 99     {
100         // 首先,自由變元有var - k個,即不確定的變元至少有var - k個.
101         for (i = k - 1; i >= 0; i--)
102         {
103             // 第i行一定不會是(0, 0, ..., 0)的情況,因為這樣的行是在第k行到第equ行.
104             // 同樣,第i行一定不會是(0, 0, ..., a), a != 0的情況,這樣的無解的.
105             free_x_num = 0; // 用於判斷該行中的不確定的變元的個數,如果超過1個,則無法求解,它們仍然為不確定的變元.
106             for (j = 0; j < var; j++)
107             {
108                 if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
109             }
110             if (free_x_num > 1) continue; // 無法求解出確定的變元.
111             // 說明就只有一個不確定的變元free_index,那么可以求解出該變元,且該變元是確定的.
112             temp = a[i][var];
113             for (j = 0; j < var; j++)
114             {
115                 if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
116             }
117             x[free_index] = temp / a[i][free_index]; // 求出該變元.
118             free_x[free_index] = 0; // 該變元是確定的.
119         }
120         return var - k; // 自由變元有var - k個.
121     }
122     // 3. 唯一解的情況: 在var * (var + 1)的增廣陣中形成嚴格的上三角陣.
123     // 計算出Xn-1, Xn-2 ... X0.
124     for (i = var - 1; i >= 0; i--)
125     {
126         temp = a[i][var];
127         for (j = i + 1; j < var; j++)
128         {
129             if (a[i][j] != 0) temp -= a[i][j] * x[j];
130         }
131         if (temp % a[i][i] != 0) return -2; // 說明有浮點數解,但無整數解.
132         x[i] = temp / a[i][i];
133     }
134 return 0;
135 }
136 
137 int main(void)
138 {
139     freopen("Input.txt", "r", stdin);
140     int i, j;
141     while (scanf("%d %d", &equ, &var) != EOF)
142     {
143         memset(a, 0, sizeof(a));
144    memset(x, 0, sizeof(x));
145    memset(free_x, 1, sizeof(free_x)); // 一開始全是不確定的變元.
146         for (i = 0; i < equ; i++)
147         {
148             for (j = 0; j < var + 1; j++)
149             {
150                 scanf("%d", &a[i][j]);
151             }
152         }
153 //        Debug();
154         free_num = Gauss();
155         if (free_num == -1) printf("無解!\n");
156    else if (free_num == -2) printf("有浮點數解,無整數解!\n");
157         else if (free_num > 0)
158         {
159             printf("無窮多解! 自由變元個數為%d\n", free_num);
160             for (i = 0; i < var; i++)
161             {
162                 if (free_x[i]) printf("x%d 是不確定的\n", i + 1);
163                 else printf("x%d: %d\n", i + 1, x[i]);
164             }
165         }
166         else
167         {
168             for (i = 0; i < var; i++)
169             {
170                 printf("x%d: %d\n", i + 1, x[i]);
171             }
172         }
173         printf("\n");
174     }
175     return 0;
176 }

參考鏈接:

[1]http://www.cnblogs.com/pegasus/archive/2011/07/31/2123195.html

[2]http://blog.sina.com.cn/s/blog_684d52a9010158ka.html

[3]http://blog.csdn.net/duanxian0621/article/details/7408887

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM