當(dāng)前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

计算矩阵的逆和行列式的值(高斯消元+LU分解)

發(fā)布時(shí)間：2023/11/30 编程问答 48 豆豆

生活随笔收集整理的這篇文章主要介紹了计算矩阵的逆和行列式的值(高斯消元+LU分解) 小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

計(jì)算矩陣的逆

選主元的高斯消元法

樸素的高斯消元法是將矩陣A和單位矩陣放在一起，通過行操作（或者列操作）將A變?yōu)閱挝痪仃?#xff0c;這個(gè)時(shí)候單位矩陣就是矩陣A的逆矩陣。從上到下將A變?yōu)樯先蔷仃嚨膹?fù)雜度為O( $n^3$ )，再?gòu)南峦蠈⑸先蔷仃囎兓癁閱挝痪仃噺?fù)雜度為O( $n^3$ )，因此總共的復(fù)雜度為O( $n^3$ ) 。

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

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

實(shí)現(xiàn)代碼：

#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){//從上往下將矩陣轉(zhuǎn)化為上三角矩陣int mark=i;for(int j=i+1;j<n;++j){//查找當(dāng)前列中最大的元素if(abs(A[mark][i]) < abs(A[j][i])){mark=j;}}if(mark != i){//如果最大元素不是當(dāng)前元素for(int j=i;j<n;++j){//對(duì)兩行元素進(jìn)行交換swap(A[i][j],A[mark][j]);}for(int j=0;j<n;++j){//對(duì)B矩陣進(jìn)行同樣的操作swap(B[i][j],B[mark][j]);}}for(int j=i+1;j<n;++j){//將后面行的第i列元素全部消去double tmp=A[j][i]/A[i][i]; //避免重復(fù)計(jì)算除數(shù)for(int k=i;k<n;++k){//A矩陣第i列前面都是０,不需要操作A[j][k]-=A[i][k]*tmp;}for(int k=0;k<n;++k){//B矩陣進(jìn)行一模一樣的操作B[j][k]-=B[i][k]*tmp;}}}for(int i=n-1;i>0;--i){//從后往前將上三角矩陣轉(zhuǎn)換為單位矩陣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對(duì)角線元素變?yōu)?#xff11;,對(duì)B進(jìn)行同樣的操作B[i][k]/=A[i][i];}}for(int k=0;k<n;++k){//對(duì)B的第一行進(jìn)行操作B[0][k]/=A[0][0];} }int main() {Matrix A,B;int n; double tmp;//讀入矩陣printf("請(qǐng)輸入矩陣的大小：");scanf("%d",&n);printf("請(qǐng)輸入%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; }

運(yùn)行結(jié)果：

LU分解法

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

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

實(shí)現(xiàn)代碼:

#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){//對(duì)A矩陣進(jìn)行LU分解for(int j=i+1;j<n;++j){//Gaussian消元A[j][i]/=A[i][i]; //將初等矩陣的值直接存放在當(dāng)前行首（因?yàn)闀?huì)變成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進(jìn)行轉(zhuǎn)置,因?yàn)槲覀兦蟮氖荁的每一列的值,我們卻是以行的方式加入數(shù)組的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("請(qǐng)輸入矩陣的大小：");scanf("%d",&n);printf("請(qǐng)輸入%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; }

運(yùn)行結(jié)果

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

計(jì)算矩陣行列式的值

采用高斯消元，將其轉(zhuǎn)換為上三角后利用行列式的性質(zhì)：上三角行列式的值等于對(duì)角線元素的乘積求出行列式的值。

實(shí)現(xiàn)代碼

#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){//從上往下將矩陣轉(zhuǎn)化為上三角矩陣int mark=i;for(int j=i+1;j<n;++j){//查找當(dāng)前列中最大的元素if(abs(A[mark][i]) < abs(A[j][i])){mark=j;}}if(mark != i){//如果最大元素不是當(dāng)前元素for(int j=i;j<n;++j) {//對(duì)兩行元素進(jìn)行交換swap(A[i][j], A[mark][j]);}}for(int j=i+1;j<n;++j){//將后面行的第i列元素全部消去double tmp=A[j][i]/A[i][i]; //避免重復(fù)計(jì)算除數(shù)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("請(qǐng)輸入行列式的大小：");scanf("%d",&n);printf("請(qǐng)輸入%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; }

運(yùn)行結(jié)果：

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

LU分解復(fù)雜度分析

單純LU分解的時(shí)間復(fù)雜度其實(shí)就是高斯消元的時(shí)間復(fù)雜度。 $T(n)=2[(n?1)n+(n?2)?(n?1)+....]=2[∑i=1ni(i?1)=∑i=1ni2?∑i=1ni]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]$ 根據(jù)求和公式 $T(n)=2[n(n+1)(2n+1)6?n(n+1)2]=2n33?2n3=O(2n33)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分解解多個(gè)系數(shù)相同的n元方程組的時(shí)候可以將復(fù)雜度均攤。

總結(jié)

以上是生活随笔為你收集整理的计算矩阵的逆和行列式的值(高斯消元+LU分解)的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。

上一篇：为什么我lol选择英雄时候有声音进去就没
下一篇：证明AVL树的上界和下界