矩阵求逆的几种方法总结(C++)
2017-02-20 20:00
507 查看
矩阵求逆运算有多种算法:
伴随矩阵的思想,分别算出其伴随矩阵和行列式,再算出逆矩阵;
LU分解法(若选主元即为LUP分解法: Ax = b ==> PAx = Pb ==>LUx = Pb ==> Ly = Pb ==> Ux = y ,每步重新选主元),它有两种不同的实现;
A-1=(LU)-1=U-1L-1,将A分解为LU后,对L和U分别求逆,再相乘;
通过解线程方程组Ax=b的方式求逆矩阵。b分别取单位阵的各个列向量,所得到的解向量x就是逆矩阵的各个列向量,拼成逆矩阵即可。
下面是这两种方法的c++代码实现,所有代码均利用常规数据集验证过。
文内程序旨在实现求逆运算核心思想,某些异常检测的功能就未实现(如矩阵维数检测、矩阵奇异等)。
注意:文中A阵均为方阵。
伴随矩阵法C++程序:
LU分解法C++程序:
两种方法运行时间测试样例(运行环境不同可能会有不同结果,我的主频是2.6GHz,内存8g。时间单位:毫秒ms)
个人认为LU分解法的两个1ms其实是不准确的(实际应远小于1ms,有兴趣可以试试看)。
三种方法的复杂度分析:
伴随矩阵法:此法的时间复杂度主要来源于计算行列式,由于计算行列式的函数为递归形式,其复杂度为O(n2)[参见这里],而整体算法需要计算每个元素的代数余子式,时间复杂度直接扩大n2倍,变为O(n4)。而递归算法本身是需要占用栈空间的,因此需要注意:当矩阵的维数较大时,随着递归深度的加大,临时占用的空间将会越来越多,甚至可能会出现栈不够用的情况(当然本次实现没有遇到,因为此时的时间开销实在令人难以忍受)!
LU分解法:此法主要是分解过程耗时,求解三角矩阵的时间复杂度是O(n2),分解过程是O(n3),总体来说和高斯消元法差不多,但是避免了高斯消元法的主元素为0的过程。 为了节省空间,A=LU分解的元素存放在A的矩阵中(因为当用过了a[i][j]元素后,便不再用了,所以可以占用原矩阵A的空间)。但是有利就有弊,考虑如果是上千个元素的矩阵,引用传参,这样就改变原矩阵了,因此程序中使用A_mintor作为副本进行使用。另外,可以看出,当矩阵维数超过某值时,内存空间便不够用了(具体是多少没有试验)。还需注意的一点是:程序中未对矩阵是否奇异进行检查,如果矩阵奇异,就不应再进行下去了。
LU分解法中,还可以先分别求出U和L的逆,再相乘,此法其实与常规LU分解法差不多。
其他:
文章中用到了矩阵的原地转置算法,具体请参考第4篇文献,这种方法降低了空间复杂度。
需要注意的问题:
本文介绍的方法new了一些指针,未释放,会出现内存泄漏,使用前请释放掉。
本文参考了以下几篇文章:
http://www.cnblogs.com/tianya2543/p/3999118.html
http://blog.csdn.net/cumtwyc/article/details/49980063
http://blog.csdn.net/xx_123_1_rj/article/details/39553809
http://www.jb51.net/article/53715.htm
伴随矩阵的思想,分别算出其伴随矩阵和行列式,再算出逆矩阵;
LU分解法(若选主元即为LUP分解法: Ax = b ==> PAx = Pb ==>LUx = Pb ==> Ly = Pb ==> Ux = y ,每步重新选主元),它有两种不同的实现;
A-1=(LU)-1=U-1L-1,将A分解为LU后,对L和U分别求逆,再相乘;
通过解线程方程组Ax=b的方式求逆矩阵。b分别取单位阵的各个列向量,所得到的解向量x就是逆矩阵的各个列向量,拼成逆矩阵即可。
下面是这两种方法的c++代码实现,所有代码均利用常规数据集验证过。
文内程序旨在实现求逆运算核心思想,某些异常检测的功能就未实现(如矩阵维数检测、矩阵奇异等)。
注意:文中A阵均为方阵。
伴随矩阵法C++程序:
#include <iostream> #include <ctime> //用于产生随机数据的种子 #define N 3 //测试矩阵维数定义 //按第一行展开计算|A| double getA(double arcs ,int n) { if(n==1) { return arcs[0][0]; } double ans = 0; double temp ={0.0}; int i,j,k; for(i=0;i<n;i++) { for(j=0;j<n-1;j++) { for(k=0;k<n-1;k++) { temp[j][k] = arcs[j+1][(k>=i)?k+1:k]; } } double t = getA(temp,n-1); if(i%2==0) { ans += arcs[0][i]*t; } else { ans -= arcs[0][i]*t; } } return ans; } //计算每一行每一列的每个元素所对应的余子式,组成A* void getAStart(double arcs ,int n,double ans ) { if(n==1) { ans[0][0] = 1; return; } int i,j,k,t; double temp ; for(i=0;i<n;i++) { for(j=0;j<n;j++) { for(k=0;k<n-1;k++) { for(t=0;t<n-1;t++) { temp[k][t] = arcs[k>=i?k+1:k][t>=j?t+1:t]; } } ans[j][i] = getA(temp,n-1); //此处顺便进行了转置 if((i+j)%2 == 1) { ans[j][i] = - ans[j][i]; } } } } //得到给定矩阵src的逆矩阵保存到des中。 bool GetMatrixInverse(double src ,int n,double des ) { double flag=getA(src,n); double t ; if(0==flag) { cout<< "原矩阵行列式为0,无法求逆。请重新运行" <<endl; return false;//如果算出矩阵的行列式为0,则不往下进行 } else { getAStart(src,n,t); for(int i=0;i<n;i++) { for(int j=0;j<n;j++) { des[i][j]=t[i][j]/flag; } } } return true; } int main() { bool flag;//标志位,如果行列式为0,则结束程序 int row =N; int col=N; double matrix_before {};//{1,2,3,4,5,6,7,8,9}; //随机数据,可替换 srand((unsigned)time(0)); for(int i=0; i<N ;i++) { for(int j=0; j<N;j++) { matrix_before[i][j]=rand()%100 *0.01; } } cout<<"原矩阵:"<<endl; for(int i=0; i<N ;i++) { for(int j=0; j<N;j++) { //cout << matrix_before[i][j] <<" "; cout << *(*(matrix_before+i)+j)<<" "; } cout<<endl; } double matrix_after {}; flag=GetMatrixInverse(matrix_before,N,matrix_after); if(false==flag) return 0; cout<<"逆矩阵:"<<endl; for(int i=0; i<row ;i++) { for(int j=0; j<col;j++) { cout <<matrix_after[i][j] <<" "; //cout << *(*(matrix_after+i)+j)<<" "; } cout<<endl; } GetMatrixInverse(matrix_after,N,matrix_before); cout<<"反算的原矩阵:"<<endl;//为了验证程序的精度 for(int i=0; i<N ;i++) { for(int j=0; j<N;j++) { //cout << matrix_before[i][j] <<" "; cout << *(*(matrix_before+i)+j)<<" "; } cout<<endl; } return 0; }
LU分解法C++程序:
#include <iostream> #include <cmath> #include <ctime> #define N 300 //矩阵乘法 double * mul(double A[N*N],double B[N*N]) { double *C=new double[N*N]{}; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { for(int k=0;k<N;k++) { C[i*N+j] += A[i*N+k]*B[k*N+j]; } } } //若绝对值小于10^-10,则置为0(这是我自己定的) for(int i=0;i<N*N;i++) { if(abs(C[i])<pow(10,-10)) { C[i]=0; } } return C; } //LUP分解 void LUP_Descomposition(double A[N*N],double L[N*N],double U[N*N],int P ) { int row=0; for(int i=0;i<N;i++) { P[i]=i; } for(int i=0;i<N-1;i++) { double p=0.0d; for(int j=i;j<N;j++) { if(abs(A[j*N+i])>p) { p=abs(A[j*N+i]); row=j; } } if(0==p) { cout<< "矩阵奇异,无法计算逆" <<endl; return ; } //交换P[i]和P[row] int tmp=P[i]; P[i]=P[row]; P[row]=tmp; double tmp2=0.0d; for(int j=0;j<N;j++) { //交换A[i][j]和 A[row][j] tmp2=A[i*N+j]; A[i*N+j]=A[row*N+j]; A[row*N+j]=tmp2; } //以下同LU分解 double u=A[i*N+i],l=0.0d; for(int j=i+1;j<N;j++) { l=A[j*N+i]/u; A[j*N+i]=l; for(int k=i+1;k<N;k++) { A[j*N+k]=A[j*N+k]-A[i*N+k]*l; } } } //构造L和U for(int i=0;i<N;i++) { for(int j=0;j<=i;j++) { if(i!=j) { L[i*N+j]=A[i*N+j]; } else { L[i*N+j]=1; } } for(int k=i;k<N;k++) { U[i*N+k]=A[i*N+k]; } } } //LUP求解方程 double * LUP_Solve(double L[N*N],double U[N*N],int P ,double b ) { double *x=new double (); double *y=new double (); //正向替换 for(int i = 0;i < N;i++) { y[i] = b[P[i]]; for(int j = 0;j < i;j++) { y[i] = y[i] - L[i*N+j]*y[j]; } } //反向替换 for(int i = N-1;i >= 0; i--) { x[i]=y[i]; for(int j = N-1;j > i;j--) { x[i] = x[i] - U[i*N+j]*x[j]; } x[i] /= U[i*N+i]; } return x; } /*****************矩阵原地转置BEGIN********************/ /* 后继 */ int getNext(int i, int m, int n) { return (i%n)*m + i/n; } /* 前驱 */ int getPre(int i, int m, int n) { return (i%m)*n + i/m; } /* 处理以下标i为起点的环 */ void movedata(double *mtx, int i, int m, int n) { double temp = mtx[i]; // 暂存 int cur = i; // 当前下标 int pre = getPre(cur, m, n); while(pre != i) { mtx[cur] = mtx[pre]; cur = pre; pre = getPre(cur, m, n); } mtx[cur] = temp; } /* 转置,即循环处理所有环 */ void transpose(double *mtx, int m, int n) { for(int i=0; i<m*n; ++i) { int next = getNext(i, m, n); while(next > i) // 若存在后继小于i说明重复,就不进行下去了(只有不重复时进入while循环) next = getNext(next, m, n); if(next == i) // 处理当前环 movedata(mtx, i, m, n); } } /*****************矩阵原地转置END********************/ //LUP求逆(将每列b求出的各列x进行组装) double * LUP_solve_inverse(double A[N*N]) { //创建矩阵A的副本,注意不能直接用A计算,因为LUP分解算法已将其改变 double *A_mirror = new double[N*N](); double *inv_A=new double[N*N]();//最终的逆矩阵(还需要转置) double *inv_A_each=new double ();//矩阵逆的各列 //double *B =new double[N*N](); double *b =new double ();//b阵为B阵的列矩阵分量 for(int i=0;i<N;i++) { double *L=new double[N*N](); double *U=new double[N*N](); int *P=new int (); //构造单位阵的每一列 for(int i=0;i<N;i++) { b[i]=0; } b[i]=1; //每次都需要重新将A复制一份 for(int i=0;i<N*N;i++) { A_mirror[i]=A[i]; } LUP_Descomposition(A_mirror,L,U,P); inv_A_each=LUP_Solve (L,U,P,b); memcpy(inv_A+i*N,inv_A_each,N*sizeof(double));//将各列拼接起来 } transpose(inv_A,N,N);//由于现在根据每列b算出的x按行存储,因此需转置 return inv_A; } int main() { double *A = new double[N*N](); srand((unsigned)time(0)); for(int i=0; i<N ;i++) { for(int j=0; j<N;j++) { A[i*N+j]=rand()%100 *0.01; } } double *E_test = new double[N*N](); double *invOfA = new double[N*N](); invOfA=LUP_solve_inverse(A); E_test=mul(A,invOfA); //验证精确度 cout<< "矩阵A:" <<endl; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { cout<< A[i*N+j]<< " " ; } cout<<endl; } cout<< "inv_A:" <<endl; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { cout<< invOfA[i*N+j]<< " " ; } cout<<endl; } cout<< "E_test:" <<endl; for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { cout<< E_test[i*N+j]<< " " ; } cout<<endl; } return 0; }
两种方法运行时间测试样例(运行环境不同可能会有不同结果,我的主频是2.6GHz,内存8g。时间单位:毫秒ms)
个人认为LU分解法的两个1ms其实是不准确的(实际应远小于1ms,有兴趣可以试试看)。
三种方法的复杂度分析:
伴随矩阵法:此法的时间复杂度主要来源于计算行列式,由于计算行列式的函数为递归形式,其复杂度为O(n2)[参见这里],而整体算法需要计算每个元素的代数余子式,时间复杂度直接扩大n2倍,变为O(n4)。而递归算法本身是需要占用栈空间的,因此需要注意:当矩阵的维数较大时,随着递归深度的加大,临时占用的空间将会越来越多,甚至可能会出现栈不够用的情况(当然本次实现没有遇到,因为此时的时间开销实在令人难以忍受)!
LU分解法:此法主要是分解过程耗时,求解三角矩阵的时间复杂度是O(n2),分解过程是O(n3),总体来说和高斯消元法差不多,但是避免了高斯消元法的主元素为0的过程。 为了节省空间,A=LU分解的元素存放在A的矩阵中(因为当用过了a[i][j]元素后,便不再用了,所以可以占用原矩阵A的空间)。但是有利就有弊,考虑如果是上千个元素的矩阵,引用传参,这样就改变原矩阵了,因此程序中使用A_mintor作为副本进行使用。另外,可以看出,当矩阵维数超过某值时,内存空间便不够用了(具体是多少没有试验)。还需注意的一点是:程序中未对矩阵是否奇异进行检查,如果矩阵奇异,就不应再进行下去了。
LU分解法中,还可以先分别求出U和L的逆,再相乘,此法其实与常规LU分解法差不多。
其他:
文章中用到了矩阵的原地转置算法,具体请参考第4篇文献,这种方法降低了空间复杂度。
需要注意的问题:
本文介绍的方法new了一些指针,未释放,会出现内存泄漏,使用前请释放掉。
本文参考了以下几篇文章:
http://www.cnblogs.com/tianya2543/p/3999118.html
http://blog.csdn.net/cumtwyc/article/details/49980063
http://blog.csdn.net/xx_123_1_rj/article/details/39553809
http://www.jb51.net/article/53715.htm
相关文章推荐
- iOS ViewController跳转界面的几种方法简单总结
- 遍历python字典几种方法总结(推荐)
- jquery插入节点的几种方法总结
- 基于并发服务器几种实现方法(总结)
- javascript中创建对象的几种方法总结
- C#执行Javascript代码的几种方法总结
- Android中Activity之间数据传递的几种方法总结
- c/c++在windows下获取时间和计算时间差的几种方法总结
- java中Object转String的几种方法总结
- java 获取字节码文件的几种方法总结
- 总结几种常见的Word转换PDF方法
- 关于iPhone/iPad全屏截图与区域截图的几种方法总结
- 总结C#中窗体间传递数据的几种方法 (由别人的方法整理)
- ArcGIS Engine中删除要素的几种方法总结
- 判断素数的几种方法的总结
- 总结几种结构体初始化方法
- 异常重抛的几种方法总结
- 几种图的构建方法的总结
- C#遍历DataSet中数据的几种方法总结