网站开发产生的材料,广州前20跨境电商公司,网页设计公司取名,响应式网站 向下兼容计算矩阵的逆 选主元的高斯消元法 朴素的高斯消元法是将矩阵A和单位矩阵放在一起#xff0c;通过行操作#xff08;或者列操作#xff09;将A变为单位矩阵#xff0c;这个时候单位矩阵就是矩阵A的逆矩阵。从上到下将A变为上三角矩阵的复杂度为O(n3n^3n3)#xff0c;再从下…计算矩阵的逆 选主元的高斯消元法 朴素的高斯消元法是将矩阵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 vectorusing namespace std;typedef vectorvectordouble Matrix;
typedef vectordouble Line;void Multi(Matrix A,Matrix B)
{int nA.size();for(int i0;in;i){for(int j0;jn;j){double tmp0;for(int k0;kn;k){tmpA[i][k]*B[k][j];}printf(%6.3f ,tmp);}printf(\n);}
}void Show(Matrix A)
{int nA.size();for(int i0;in;i){for(int j0;jn;j){printf(%6.3f ,A[i][j]);}printf(\n);}
}void Gaussian(Matrix A,Matrix B)
{int nA.size();Line x(n,0); x[0]1; B.push_back(x);for(int i1;in;i){//将B初始化为单位矩阵x[i-1]0; x[i]1;B.push_back(x);}for(int i0;in-1;i){//从上往下将矩阵转化为上三角矩阵int marki;for(int ji1;jn;j){//查找当前列中最大的元素if(abs(A[mark][i]) abs(A[j][i])){markj;}}if(mark ! i){//如果最大元素不是当前元素for(int ji;jn;j){//对两行元素进行交换swap(A[i][j],A[mark][j]);}for(int j0;jn;j){//对B矩阵进行同样的操作swap(B[i][j],B[mark][j]);}}for(int ji1;jn;j){//将后面行的第i列元素全部消去double tmpA[j][i]/A[i][i]; //避免重复计算除数for(int ki;kn;k){//A矩阵第i列前面都是,不需要操作A[j][k]-A[i][k]*tmp;}for(int k0;kn;k){//B矩阵进行一模一样的操作B[j][k]-B[i][k]*tmp;}}}for(int in-1;i0;--i){//从后往前将上三角矩阵转换为单位矩阵for(int ji-1;j0;--j){//将前面行的第i列元素全部消去double tmpA[j][i]/A[i][i];A[j][i]0;//其他列的元素不变for(int kn-1;k0;--k){B[j][k]-B[i][k]*tmp;}}for(int k0;kn;k){//将A对角线元素变为,对B进行同样的操作B[i][k]/A[i][i];}}for(int k0;kn;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 i0;in;i){Line x;A.push_back(x);for(int j0;jn;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(ij)Ri-Rj*C (ij)Ri−Rj∗C(ij)相当于左乘初等矩阵P[i,j]−cP[i,j]-cP[i,j]−c我们将这些信息保存下来就能得到Pn−1∗Pn−2∗...∗P1∗UAP_{n-1}*P_{n-2}*...*P_1*UAPn−1∗Pn−2∗...∗P1∗UA其中U是上三角矩阵也就是我们高斯消元以后得到的矩阵。根据矩阵运算的结合律我们将初等矩阵合并为矩阵LLL得到L∗UAL*UAL∗UA这里面的L,UL,UL,U均为三角矩阵然后利用三角矩阵解方程组将会非常容易。例如我们需要求解AXBAXBAXB即就是求解LUXBLUXBLUXB我们令YUXYUXYUX则解两个含有三角阵的方程就可以解决问题。 观察L,UL,UL,U矩阵我们发现可以将他们合并而且可以在原矩阵中原地操作因此空间复杂度为O(1)O(1)O(1)。而且对矩阵的LU分解可以重复多次运算我们可以借此将对矩阵求逆的过程转换为AA−1EAA^{-1}EAA−1E求解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 vectorusing namespace std;typedef vectorvectordouble Matrix;
typedef vectordouble Line;void Multi(Matrix A,Matrix B)
{int nA.size();for(int i0;in;i){for(int j0;jn;j){double tmp0;for(int k0;kn;k){tmpA[i][k]*B[k][j];}printf(%6.3f ,tmp);}printf(\n);}
}void Show(Matrix A)
{int nA.size();for(int i0;in;i){for(int j0;jn;j){printf(%6.3f ,A[i][j]);}printf(\n);}
}Line LUWork(Matrix A,Line Z)
{int nA.size();//解LYZLine Y(n);for(int i0;in;i){double sum0;for(int j0;ji;j){sumA[i][j]*Y[j];}Y[i]Z[i]-sum;}//解UXYLine X(n);for(int in-1;i0;--i){double sum0;for(int jn-1;ji;--j){sumA[i][j]*X[j];}X[i](Y[i]-sum)/A[i][i];}return X;
}void LU(Matrix A,Matrix B)
{int nA.size();for(int i0;in-1;i){//对A矩阵进行LU分解for(int ji1;jn;j){//Gaussian消元A[j][i]/A[i][i]; //将初等矩阵的值直接存放在当前行首因为会变成0,没有什么作用for(int ki1;kn;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 i1;in;i){Z[i-1]0; Z[i]1;B.push_back(LUWork(A,Z));}//将B进行转置,因为我们求的是B的每一列的值,我们却是以行的方式加入数组的for(int i0;in;i){for(int j0;ji;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 i0;in;i){Line x;A.push_back(x);for(int j0;jn;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−1A∗∣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 vectorusing namespace std;typedef vectorvectordouble Matrix;
typedef vectordouble Line;double Gaussian(Matrix A)
{int nA.size();for(int i0;in-1;i){//从上往下将矩阵转化为上三角矩阵int marki;for(int ji1;jn;j){//查找当前列中最大的元素if(abs(A[mark][i]) abs(A[j][i])){markj;}}if(mark ! i){//如果最大元素不是当前元素for(int ji;jn;j) {//对两行元素进行交换swap(A[i][j], A[mark][j]);}}for(int ji1;jn;j){//将后面行的第i列元素全部消去double tmpA[j][i]/A[i][i]; //避免重复计算除数for(int ki;kn;k){//A矩阵第i列前面都是,不需要操作A[j][k]-A[i][k]*tmp;}}}double ret1.0;for(int i0;in;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 i0;in;i){Line x;A.push_back(x);for(int j0;jn;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[∑i1ni(i−1)∑i1ni2−∑i1ni]T(n)2[(n-1)n(n-2)*(n-1)....]2[\sum_{i1}^n i(i-1)\sum_{i1}^n i^2 - \sum_{i1}^n i]T(n)2[(n−1)n(n−2)∗(n−1)....]2[∑i1ni(i−1)∑i1ni2−∑i1ni]根据求和公式T(n)2[n(n1)(2n1)6−n(n1)2]2n33−2n3O(2n33)T(n)2[\frac{n(n1)(2n1)}{6}-\frac{n(n1)}{2}]\frac{2n^3}{3}-\frac{2n}{3}O(\frac{2n^3}{3})T(n)2[6n(n1)(2n1)−2n(n1)]32n3−32nO(32n3)。
如果使用LU分解解多个系数相同的n元方程组的时候可以将复杂度均摊。
