計算矩陣的逆

1. 選主元的高斯消元法

   樸素的高斯消元法是將矩陣A和單位矩陣放在一起，通過行操作（或者列操作）將A變為單位矩陣，這個時候單位矩陣就是矩陣A的逆矩陣。從上到下將A變為上三角矩陣的復雜度為O($n^3$)，再從下往上將上三角矩陣變化為單位矩陣復雜度為O($n^3$)，因此總共的復雜度為O($n^3$) 。

   還有一種做法是按照高斯消元接線性方程組的方式求解n次線性方程組，這樣復雜度為O($n^4$)，因此我們采用上面的方法。

   上面的方法雖然可以有效的解決問題，但是在計算機中計算除法的時候如果除數太小將會產生較大的誤差，因此我們就需要主動采取措施。一個簡單的方法就是選擇主元，即找一個最大的每次將要去消除其他行的行首元素。根據查找范圍的不同又分為全選主元和列選主元兩種。全選主元顧名思義就是在當前行下方所有元素中尋找最大的元素，然后將其行和列置換到合適的位置。這個操作是O($n^2$)的，但是可以有效避免除數過小的問題。這里簡單起見采用的是列選主元：在當前列中選擇一個最大的元素將其置換到當前行。

   實現代碼：

   #include <iostream>
   #include <algorithm>
   #include <vector>using namespace std;typedef vector<vector<double> > Matrix;
   typedef vector<double> Line;void Multi(Matrix &A,Matrix &B)
   {int n=A.size();for(int i=0;i<n;++i){for(int j=0;j<n;++j){double tmp=0;for(int k=0;k<n;++k){tmp+=A[i][k]*B[k][j];}printf("%6.3f ",tmp);}printf("\n");}
   }void Show(Matrix &A)
   {int n=A.size();for(int i=0;i<n;++i){for(int j=0;j<n;++j){printf("%6.3f ",A[i][j]);}printf("\n");}
   }void Gaussian(Matrix A,Matrix& B)
   {int n=A.size();Line x(n,0); x[0]=1;  B.push_back(x);for(int i=1;i<n;++i){//將B初始化為單位矩陣x[i-1]=0; x[i]=1;B.push_back(x);}for(int i=0;i<n-1;++i){//從上往下將矩陣轉化為上三角矩陣int mark=i;for(int j=i+1;j<n;++j){//查找當前列中最大的元素if(abs(A[mark][i]) < abs(A[j][i])){mark=j;}}if(mark != i){//如果最大元素不是當前元素for(int j=i;j<n;++j){//對兩行元素進行交換swap(A[i][j],A[mark][j]);}for(int j=0;j<n;++j){//對B矩陣進行同樣的操作swap(B[i][j],B[mark][j]);}}for(int j=i+1;j<n;++j){//將后面行的第i列元素全部消去double tmp=A[j][i]/A[i][i]; //避免重復計算除數for(int k=i;k<n;++k){//A矩陣第i列前面都是０,不需要操作A[j][k]-=A[i][k]*tmp;}for(int k=0;k<n;++k){//B矩陣進行一模一樣的操作B[j][k]-=B[i][k]*tmp;}}}for(int i=n-1;i>0;--i){//從后往前將上三角矩陣轉換為單位矩陣for(int j=i-1;j>=0;--j){//將前面行的第i列元素全部消去double tmp=A[j][i]/A[i][i];A[j][i]=0;//其他列的元素不變for(int k=n-1;k>=0;--k){B[j][k]-=B[i][k]*tmp;}}for(int k=0;k<n;++k){//將A對角線元素變為１,對B進行同樣的操作B[i][k]/=A[i][i];}}for(int k=0;k<n;++k){//對B的第一行進行操作B[0][k]/=A[0][0];}
   }int main()
   {Matrix A,B;int n; double tmp;//讀入矩陣printf("請輸入矩陣的大小：");scanf("%d",&n);printf("請輸入%d*%d的矩陣：\n",n,n);for(int i=0;i<n;++i){Line x;A.push_back(x);for(int j=0;j<n;++j){scanf("%lf",&tmp);A[i].push_back(tmp);}}Gaussian(A,B);printf("Inverse Matrix:\n");Show(B);printf("A X B:\n");Multi(A,B);return 0;
   }
   

   運行結果：

   

2. LU分解法

   LU分解可以看做是高斯消元的一個變化的應用，不同的地方在于它將每次高斯消元的步驟都保存了下來。可以這樣做的原因是我們對矩陣的行操作都可以看做我們給矩陣左乘了一個矩陣。例如，我們在高斯消元的過程中有如下操作：$Ri?Rj?C(i>j)$，相當于左乘初等矩陣$P[i,j]=?c$，我們將這些信息保存下來就能得到$P_{n?1}?P_{n?2}?...?P_1?U=A$，其中U是上三角矩陣，也就是我們高斯消元以后得到的矩陣。根據矩陣運算的結合律我們將初等矩陣合并為矩陣$L$，得到$L?U=A$，這里面的$L,U$均為三角矩陣，然后利用三角矩陣解方程組將會非常容易。例如我們需要求解$AX=B$，即就是求解$LUX=B$，我們令$Y=UX$，則解兩個含有三角陣的方程就可以解決問題。

   觀察$L,U$矩陣，我們發現可以將他們合并，而且可以在原矩陣中原地操作，因此空間復雜度為$O(1)$。而且對矩陣的LU分解可以重復多次運算，我們可以借此將對矩陣求逆的過程轉換為$A{A}^{?1}=E$，求解A的n個n元方程組。復雜度為LU分解$O(n^3)$加上n次求解方程組$n?O(n^2)$，因此總共的復雜度為$O(n^3)$。

   實現代碼:

   #include <iostream>
   #include <algorithm>
   #include <vector>using namespace std;typedef vector<vector<double> > Matrix;
   typedef vector<double> Line;void Multi(Matrix &A,Matrix &B)
   {int n=A.size();for(int i=0;i<n;++i){for(int j=0;j<n;++j){double tmp=0;for(int k=0;k<n;++k){tmp+=A[i][k]*B[k][j];}printf("%6.3f ",tmp);}printf("\n");}
   }void Show(Matrix &A)
   {int n=A.size();for(int i=0;i<n;++i){for(int j=0;j<n;++j){printf("%6.3f ",A[i][j]);}printf("\n");}
   }Line LUWork(Matrix& A,Line& Z)
   {int n=A.size();//解LY=ZLine Y(n);for(int i=0;i<n;++i){double sum=0;for(int j=0;j<i;++j){sum+=A[i][j]*Y[j];}Y[i]=Z[i]-sum;}//解UX=YLine X(n);for(int i=n-1;i>=0;--i){double sum=0;for(int j=n-1;j>i;--j){sum+=A[i][j]*X[j];}X[i]=(Y[i]-sum)/A[i][i];}return X;
   }void LU(Matrix A,Matrix &B)
   {int n=A.size();for(int i=0;i<n-1;++i){//對A矩陣進行LU分解for(int j=i+1;j<n;++j){//Gaussian消元A[j][i]/=A[i][i]; //將初等矩陣的值直接存放在當前行首（因為會變成0,沒有什么作用）for(int k=i+1;k<n;++k){A[j][k]-=A[i][k]*A[j][i];}}}//解n次n元方程組Line Z(n,0); Z[0]=1; B.push_back(LUWork(A,Z));for(int i=1;i<n;++i){Z[i-1]=0; Z[i]=1;B.push_back(LUWork(A,Z));}//將B進行轉置,因為我們求的是B的每一列的值,我們卻是以行的方式加入數組的for(int i=0;i<n;++i){for(int j=0;j<i;++j){swap(B[i][j],B[j][i]);}}
   }int main()
   {Matrix A,B;int n; double tmp;//讀入矩陣printf("請輸入矩陣的大小：");scanf("%d",&n);printf("請輸入%d*%d的矩陣：\n",n,n);for(int i=0;i<n;++i){Line x;A.push_back(x);for(int j=0;j<n;++j){scanf("%lf",&tmp);A[i].push_back(tmp);}}LU(A,B);printf("Inverse Matrix:\n");Show(B);printf("A X B:\n");Multi(A,B);return 0;
   }

   運行結果

   

3. 我們還可以通過$A^{?1}=\frac{A^{?}}{|A|}$來求矩陣的逆，但是通過下面的討論我們發現求|A|的復雜度至少也是$O(n^3)$的，還需要求解$A^{?}$。因此這種方式不實用。

---

計算矩陣行列式的值

1. 采用高斯消元，將其轉換為上三角后利用行列式的性質：上三角行列式的值等于對角線元素的乘積求出行列式的值。

   實現代碼

   #include <iostream>
   #include <algorithm>
   #include <vector>using namespace std;typedef vector<vector<double> > Matrix;
   typedef vector<double> Line;double Gaussian(Matrix A)
   {int n=A.size();for(int i=0;i<n-1;++i){//從上往下將矩陣轉化為上三角矩陣int mark=i;for(int j=i+1;j<n;++j){//查找當前列中最大的元素if(abs(A[mark][i]) < abs(A[j][i])){mark=j;}}if(mark != i){//如果最大元素不是當前元素for(int j=i;j<n;++j) {//對兩行元素進行交換swap(A[i][j], A[mark][j]);}}for(int j=i+1;j<n;++j){//將后面行的第i列元素全部消去double tmp=A[j][i]/A[i][i]; //避免重復計算除數for(int k=i;k<n;++k){//A矩陣第i列前面都是０,不需要操作A[j][k]-=A[i][k]*tmp;}}}double ret=1.0;for(int i=0;i<n;++i){ret*=A[i][i];}return ret;
   }int main()
   {Matrix A,B;int n; double tmp;//讀入行列式printf("請輸入行列式的大小：");scanf("%d",&n);printf("請輸入%d*%d的行列式：\n",n,n);for(int i=0;i<n;++i){Line x;A.push_back(x);for(int j=0;j<n;++j){scanf("%lf",&tmp);A[i].push_back(tmp);}}printf("行列式的值為:%.f\n",Gaussian(A));return 0;
   }

   運行結果：

   

2. 我們還可以利用行列式的性質：行列式等于任意行（列）各個元素與其代數余子式的乘積的和。這樣進行遞歸求解，得到遞歸式$T(n)=nT(n?1)+n$，復雜度為$O(n!)$

---

LU分解復雜度分析

單純LU分解的時間復雜度其實就是高斯消元的時間復雜度。$T(n)=2[(n?1)n+(n?2)?(n?1)+....]=2[{\sum }_{i=1}^{n}i(i?1)={\sum }_{i=1}^n{i}^2?{\sum }_{i=1}^{n}i]$根據求和公式$T(n)=2[\frac{n(n+1)(2n+1)}{6}?\frac{n(n+1)}{2}]=\frac{2n^3}{3}?\frac{2n}{3}=O(\frac{2n^3}{3})$。

如果使用LU分解解多個系數相同的n元方程組的時候可以將復雜度均攤。

---