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

计算矩阵的逆

选主元的高斯消元法

朴素的高斯消元法是将矩阵A和单位矩阵放在一起，通过行操作（或者列操作）将A变为单位矩阵，这个时候单位矩阵就是矩阵A的逆矩阵。从上到下将A变为上三角矩阵的复杂度为O(n3n^3n3)，再从下往上将上三角矩阵变化为单位矩阵复杂度为O(n3n^3n3)，因此总共的复杂度为O(n3n^3n3) 。

还有一种做法是按照高斯消元接线性方程组的方式求解n次线性方程组，这样复杂度为O(n4n^4n4)，因此我们采用上面的方法。

上面的方法虽然可以有效的解决问题，但是在计算机中计算除法的时候如果除数太小将会产生较大的误差，因此我们就需要主动采取措施。一个简单的方法就是选择主元，即找一个最大的每次将要去消除其他行的行首元素。根据查找范围的不同又分为全选主元和列选主元两种。全选主元顾名思义就是在当前行下方所有元素中寻找最大的元素，然后将其行和列置换到合适的位置。这个操作是O(n2n^2n2)的，但是可以有效避免除数过小的问题。这里简单起见采用的是列选主元：在当前列中选择一个最大的元素将其置换到当前行。

实现代码：

#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;
}

运行结果：

LU分解法

LU分解可以看做是高斯消元的一个变化的应用，不同的地方在于它将每次高斯消元的步骤都保存了下来。可以这样做的原因是我们对矩阵的行操作都可以看做我们给矩阵左乘了一个矩阵。例如，我们在高斯消元的过程中有如下操作：Ri−Rj∗C(i>j)Ri-Rj*C (i>j)Ri−Rj∗C(i>j)，相当于左乘初等矩阵P[i,j]=−cP[i,j]=-cP[i,j]=−c，我们将这些信息保存下来就能得到Pn−1∗Pn−2∗...∗P1∗U=AP_{n-1}*P_{n-2}*...*P_1*U=APn−1∗Pn−2∗...∗P1∗U=A，其中U是上三角矩阵，也就是我们高斯消元以后得到的矩阵。根据矩阵运算的结合律我们将初等矩阵合并为矩阵LLL，得到L∗U=AL*U=AL∗U=A，这里面的L,UL,UL,U均为三角矩阵，然后利用三角矩阵解方程组将会非常容易。例如我们需要求解AX=BAX=BAX=B，即就是求解LUX=BLUX=BLUX=B，我们令Y=UXY=UXY=UX，则解两个含有三角阵的方程就可以解决问题。

观察L,UL,UL,U矩阵，我们发现可以将他们合并，而且可以在原矩阵中原地操作，因此空间复杂度为O(1)O(1)O(1)。而且对矩阵的LU分解可以重复多次运算，我们可以借此将对矩阵求逆的过程转换为AA−1=EAA^{-1}=EAA−1=E，求解A的n个n元方程组。复杂度为LU分解O(n3)O(n^3)O(n3)加上n次求解方程组n∗O(n2)n*O(n^2)n∗O(n2)，因此总共的复杂度为O(n3)O(n^3)O(n3)。

实现代码:

#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;
}

运行结果

我们还可以通过A−1=A∗∣A∣A^{-1}=\frac{A^*}{|A|}A−1=∣A∣A∗来求矩阵的逆，但是通过下面的讨论我们发现求|A|的复杂度至少也是O(n3)O(n^3)O(n3)的，还需要求解A∗A^*A∗。因此这种方式不实用。

计算矩阵行列式的值

采用高斯消元，将其转换为上三角后利用行列式的性质：上三角行列式的值等于对角线元素的乘积求出行列式的值。

实现代码

#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;
}

运行结果：

我们还可以利用行列式的性质：行列式等于任意行（列）各个元素与其代数余子式的乘积的和。这样进行递归求解，得到递归式T(n)=nT(n−1)+nT(n)=nT(n-1)+nT(n)=nT(n−1)+n，复杂度为O(n!)O(n!)O(n!)

LU分解复杂度分析

单纯LU分解的时间复杂度其实就是高斯消元的时间复杂度。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]T(n)=2[(n−1)n+(n−2)∗(n−1)+....]=2[∑i=1ni(i−1)=∑i=1ni2−∑i=1ni]根据求和公式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})T(n)=2[6n(n+1)(2n+1)−2n(n+1)]=32n3−32n=O(32n3)。

如果使用LU分解解多个系数相同的n元方程组的时候可以将复杂度均摊。