矩阵求逆的几种方法总结(C++)
矩阵求逆运算有多种算法:
- 伴随矩阵的思想,分别算出其伴随矩阵和行列式,再算出逆矩阵;
- 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++程序:
1 #include <iostream> 2 #include <ctime> //用于产生随机数据的种子 3 4 #define N 3 //测试矩阵维数定义 5 6 //按第一行展开计算|A| 7 double getA(double arcs[N][N],int n) 8 { 9 if(n==1) 10 { 11 return arcs[0][0]; 12 } 13 double ans = 0; 14 double temp[N][N]={0.0}; 15 int i,j,k; 16 for(i=0;i<n;i++) 17 { 18 for(j=0;j<n-1;j++) 19 { 20 for(k=0;k<n-1;k++) 21 { 22 temp[j][k] = arcs[j+1][(k>=i)?k+1:k]; 23 24 } 25 } 26 double t = getA(temp,n-1); 27 if(i%2==0) 28 { 29 ans += arcs[0][i]*t; 30 } 31 else 32 { 33 ans -= arcs[0][i]*t; 34 } 35 } 36 return ans; 37 } 38 39 //计算每一行每一列的每个元素所对应的余子式,组成A* 40 void getAStart(double arcs[N][N],int n,double ans[N][N]) 41 { 42 if(n==1) 43 { 44 ans[0][0] = 1; 45 return; 46 } 47 int i,j,k,t; 48 double temp[N][N]; 49 for(i=0;i<n;i++) 50 { 51 for(j=0;j<n;j++) 52 { 53 for(k=0;k<n-1;k++) 54 { 55 for(t=0;t<n-1;t++) 56 { 57 temp[k][t] = arcs[k>=i?k+1:k][t>=j?t+1:t]; 58 } 59 } 60 61 62 ans[j][i] = getA(temp,n-1); //此处顺便进行了转置 63 if((i+j)%2 == 1) 64 { 65 ans[j][i] = - ans[j][i]; 66 } 67 } 68 } 69 } 70 71 //得到给定矩阵src的逆矩阵保存到des中。 72 bool GetMatrixInverse(double src[N][N],int n,double des[N][N]) 73 { 74 double flag=getA(src,n); 75 double t[N][N]; 76 if(0==flag) 77 { 78 cout<< "原矩阵行列式为0,无法求逆。请重新运行" <<endl; 79 return false;//如果算出矩阵的行列式为0,则不往下进行 80 } 81 else 82 { 83 getAStart(src,n,t); 84 for(int i=0;i<n;i++) 85 { 86 for(int j=0;j<n;j++) 87 { 88 des[i][j]=t[i][j]/flag; 89 } 90 91 } 92 } 93 94 return true; 95 } 96 97 int main() 98 { 99 bool flag;//标志位,如果行列式为0,则结束程序 100 int row =N; 101 int col=N; 102 double matrix_before[N][N]{};//{1,2,3,4,5,6,7,8,9}; 103 104 //随机数据,可替换 105 srand((unsigned)time(0)); 106 for(int i=0; i<N ;i++) 107 { 108 for(int j=0; j<N;j++) 109 { 110 matrix_before[i][j]=rand()%100 *0.01; 111 } 112 } 113 114 cout<<"原矩阵:"<<endl; 115 116 for(int i=0; i<N ;i++) 117 { 118 for(int j=0; j<N;j++) 119 { 120 //cout << matrix_before[i][j] <<" "; 121 cout << *(*(matrix_before+i)+j)<<" "; 122 } 123 cout<<endl; 124 } 125 126 127 double matrix_after[N][N]{}; 128 flag=GetMatrixInverse(matrix_before,N,matrix_after); 129 if(false==flag) 130 return 0; 131 132 133 cout<<"逆矩阵:"<<endl; 134 135 for(int i=0; i<row ;i++) 136 { 137 for(int j=0; j<col;j++) 138 { 139 cout <<matrix_after[i][j] <<" "; 140 //cout << *(*(matrix_after+i)+j)<<" "; 141 } 142 cout<<endl; 143 } 144 145 GetMatrixInverse(matrix_after,N,matrix_before); 146 147 cout<<"反算的原矩阵:"<<endl;//为了验证程序的精度 148 149 for(int i=0; i<N ;i++) 150 { 151 for(int j=0; j<N;j++) 152 { 153 //cout << matrix_before[i][j] <<" "; 154 cout << *(*(matrix_before+i)+j)<<" "; 155 } 156 cout<<endl; 157 } 158 159 160 return 0; 161 }
LU分解法C++程序:
1 #include <iostream> 2 #include <cmath> 3 #include <ctime> 4 5 #define N 300 6 7 //矩阵乘法 8 double * mul(double A[N*N],double B[N*N]) 9 { 10 double *C=new double[N*N]{}; 11 for(int i=0;i<N;i++) 12 { 13 for(int j=0;j<N;j++) 14 { 15 for(int k=0;k<N;k++) 16 { 17 C[i*N+j] += A[i*N+k]*B[k*N+j]; 18 } 19 } 20 } 21 22 //若绝对值小于10^-10,则置为0(这是我自己定的) 23 for(int i=0;i<N*N;i++) 24 { 25 if(abs(C[i])<pow(10,-10)) 26 { 27 C[i]=0; 28 } 29 } 30 31 return C; 32 } 33 34 //LUP分解 35 void LUP_Descomposition(double A[N*N],double L[N*N],double U[N*N],int P[N]) 36 { 37 int row=0; 38 for(int i=0;i<N;i++) 39 { 40 P[i]=i; 41 } 42 for(int i=0;i<N-1;i++) 43 { 44 double p=0.0d; 45 for(int j=i;j<N;j++) 46 { 47 if(abs(A[j*N+i])>p) 48 { 49 p=abs(A[j*N+i]); 50 row=j; 51 } 52 } 53 if(0==p) 54 { 55 cout<< "矩阵奇异,无法计算逆" <<endl; 56 return ; 57 } 58 59 //交换P[i]和P[row] 60 int tmp=P[i]; 61 P[i]=P[row]; 62 P[row]=tmp; 63 64 double tmp2=0.0d; 65 for(int j=0;j<N;j++) 66 { 67 //交换A[i][j]和 A[row][j] 68 tmp2=A[i*N+j]; 69 A[i*N+j]=A[row*N+j]; 70 A[row*N+j]=tmp2; 71 } 72 73 //以下同LU分解 74 double u=A[i*N+i],l=0.0d; 75 for(int j=i+1;j<N;j++) 76 { 77 l=A[j*N+i]/u; 78 A[j*N+i]=l; 79 for(int k=i+1;k<N;k++) 80 { 81 A[j*N+k]=A[j*N+k]-A[i*N+k]*l; 82 } 83 } 84 85 } 86 87 //构造L和U 88 for(int i=0;i<N;i++) 89 { 90 for(int j=0;j<=i;j++) 91 { 92 if(i!=j) 93 { 94 L[i*N+j]=A[i*N+j]; 95 } 96 else 97 { 98 L[i*N+j]=1; 99 } 100 } 101 for(int k=i;k<N;k++) 102 { 103 U[i*N+k]=A[i*N+k]; 104 } 105 } 106 107 } 108 109 //LUP求解方程 110 double * LUP_Solve(double L[N*N],double U[N*N],int P[N],double b[N]) 111 { 112 double *x=new double[N](); 113 double *y=new double[N](); 114 115 //正向替换 116 for(int i = 0;i < N;i++) 117 { 118 y[i] = b[P[i]]; 119 for(int j = 0;j < i;j++) 120 { 121 y[i] = y[i] - L[i*N+j]*y[j]; 122 } 123 } 124 //反向替换 125 for(int i = N-1;i >= 0; i--) 126 { 127 x[i]=y[i]; 128 for(int j = N-1;j > i;j--) 129 { 130 x[i] = x[i] - U[i*N+j]*x[j]; 131 } 132 x[i] /= U[i*N+i]; 133 } 134 return x; 135 } 136 137 /*****************矩阵原地转置BEGIN********************/ 138 139 /* 后继 */ 140 int getNext(int i, int m, int n) 141 { 142 return (i%n)*m + i/n; 143 } 144 145 /* 前驱 */ 146 int getPre(int i, int m, int n) 147 { 148 return (i%m)*n + i/m; 149 } 150 151 /* 处理以下标i为起点的环 */ 152 void movedata(double *mtx, int i, int m, int n) 153 { 154 double temp = mtx[i]; // 暂存 155 int cur = i; // 当前下标 156 int pre = getPre(cur, m, n); 157 while(pre != i) 158 { 159 mtx[cur] = mtx[pre]; 160 cur = pre; 161 pre = getPre(cur, m, n); 162 } 163 mtx[cur] = temp; 164 } 165 166 /* 转置,即循环处理所有环 */ 167 void transpose(double *mtx, int m, int n) 168 { 169 for(int i=0; i<m*n; ++i) 170 { 171 int next = getNext(i, m, n); 172 while(next > i) // 若存在后继小于i说明重复,就不进行下去了(只有不重复时进入while循环) 173 next = getNext(next, m, n); 174 if(next == i) // 处理当前环 175 movedata(mtx, i, m, n); 176 } 177 } 178 /*****************矩阵原地转置END********************/ 179 180 //LUP求逆(将每列b求出的各列x进行组装) 181 double * LUP_solve_inverse(double A[N*N]) 182 { 183 //创建矩阵A的副本,注意不能直接用A计算,因为LUP分解算法已将其改变 184 double *A_mirror = new double[N*N](); 185 double *inv_A=new double[N*N]();//最终的逆矩阵(还需要转置) 186 double *inv_A_each=new double[N]();//矩阵逆的各列 187 //double *B =new double[N*N](); 188 double *b =new double[N]();//b阵为B阵的列矩阵分量 189 190 for(int i=0;i<N;i++) 191 { 192 double *L=new double[N*N](); 193 double *U=new double[N*N](); 194 int *P=new int[N](); 195 196 //构造单位阵的每一列 197 for(int i=0;i<N;i++) 198 { 199 b[i]=0; 200 } 201 b[i]=1; 202 203 //每次都需要重新将A复制一份 204 for(int i=0;i<N*N;i++) 205 { 206 A_mirror[i]=A[i]; 207 } 208 209 LUP_Descomposition(A_mirror,L,U,P); 210 211 inv_A_each=LUP_Solve (L,U,P,b); 212 memcpy(inv_A+i*N,inv_A_each,N*sizeof(double));//将各列拼接起来 213 } 214 transpose(inv_A,N,N);//由于现在根据每列b算出的x按行存储,因此需转置 215 216 return inv_A; 217 } 218 219 int main() 220 { 221 double *A = new double[N*N](); 222 223 srand((unsigned)time(0)); 224 for(int i=0; i<N ;i++) 225 { 226 for(int j=0; j<N;j++) 227 { 228 A[i*N+j]=rand()%100 *0.01; 229 } 230 } 231 232 233 double *E_test = new double[N*N](); 234 double *invOfA = new double[N*N](); 235 invOfA=LUP_solve_inverse(A); 236 237 E_test=mul(A,invOfA); //验证精确度 238 239 cout<< "矩阵A:" <<endl; 240 for(int i=0;i<N;i++) 241 { 242 for(int j=0;j<N;j++) 243 { 244 cout<< A[i*N+j]<< " " ; 245 } 246 cout<<endl; 247 } 248 249 cout<< "inv_A:" <<endl; 250 for(int i=0;i<N;i++) 251 { 252 for(int j=0;j<N;j++) 253 { 254 cout<< invOfA[i*N+j]<< " " ; 255 } 256 cout<<endl; 257 } 258 259 cout<< "E_test:" <<endl; 260 for(int i=0;i<N;i++) 261 { 262 for(int j=0;j<N;j++) 263 { 264 cout<< E_test[i*N+j]<< " " ; 265 } 266 cout<<endl; 267 } 268 269 return 0; 270 }
两种方法运行时间测试样例(运行环境不同可能会有不同结果,我的主频是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
转载于:https://www.cnblogs.com/xiaoxi666/p/6421228.html
矩阵求逆的几种方法总结(C++)相关推荐
- 8种方法用Python实现线性回归,为你解析最高效选择
来源:大数据文摘 编译:丁慧.katherine Hou.钱天培 作者:TirthajyotiSarkar 本文共1856字,建议阅读6分钟. 本文为大家对比了8种用Python实现线性回归的方法哪个 ...
- spring boot项目 中止运行 最常用的几种方法
spring boot项目 中止运行 最常用的几种方法: 1. 调用接口,停止应用上下文 @RestController public class ShutdownController impleme ...
- 设置select下拉框不可修改的→“四”←种方法
设置select下拉框为不可修改的几种方法: 因为select的特殊性,导致它不能像input表单一样简单地设置一个readonly来限制修改,所以,我们需要进行别的操作! 1.为下拉框添加样式,可以 ...
- 用python下载文件的若干种方法汇总
压缩文件可以直接放到下载器里面下载的 you-get 连接 下载任意文件 重点 用python下载文件的若干种方法汇总 写文章 用python下载文件的若干种方法汇总 zhangqibot发表于Met ...
- Python两个字典键同值相加的几种方法
版权声明:本文为博主原创文章,遵循 CC 4.0 by-sa 版权协议,转载请附上原文出处链接和本声明. 本文链接:https://blog.csdn.net/Jerry_1126/article/d ...
- VS中C#读取app.config数据库配置字符串的三种方法(转)
关于VS2008或VS2005中数据库配置字符串的三种取法 VS2008建立Form程序时,如果添加数据源会在配置文件 app.config中自动写入连接字符串,这个字符串将会在你利用DataSet, ...
- 在PHP中使用全局变量的几种方法
简介 即使开发一个新的大型PHP程序,你也不可避免的要使用到全局数据,因为有些数据是需要用到你的代码的不同部分的.一些常见的全局数据有:程序设定类.数据库连接类.用户资料等等.有很多方法能够使这些数据 ...
- SQL Server中灾难时备份结尾日志(Tail of log)的两种方法
简介 在数据库数据文件因各种原因发生损坏时,如果日志文件没有损坏.可以通过备份结尾日志(Tail of log)使得数据库可以恢复到灾难发生时的状态. 例如: 上图中.在DB_1中做了完整备份,在Lo ...
- c#中分割字符串的几种方法
第一种方法:打开vs.net新建一个控制台项目.然后在Main()方法下输入下面的程序. string s="abcdeabcdeabcde"; string[] sArray=s ...
- java数据输入的步骤_Java学习日志1.4 Scanner 数据输入的三种方法
Scanner sc = new Scanner(System.in); /注意in 是InputStream的缩写,是字节输入流的意思. 整句话的含义就是: new 一个对象,接受从键盘输入的数据, ...
最新文章
- Winform中实现序列化指定类型的对象到指定的Xml文件和从指定的Xml文件中反序列化指定类型的对象
- NFL discussion调研
- IDC机房KVM应用案例分析
- c语言斐波那契数列前20项每行5个数,求c++:源程序。前20项斐波那契数列 ,要求输出的时候每行输出五个...
- 模拟电子技术基础简明课程(第三版)思维导图
- Javashop电商系统7.1.5源码发布
- 【工作提效】PLSQL使用技巧
- 6大最常用的Java机器学习库一览
- Echarts——自定义仪表盘图表
- 每周分享第 38 期
- t420i升级固态硬盘提升_旧电脑升级!使用固态硬盘必做的5件事,让win10操作流畅如win7...
- 我的世界服务器修改npc指令,我的世界自定义npc指令 | 手游网游页游攻略大全
- 企业资源计划-MPS计算(附详细解题步骤及计算过程)
- 查看并彻底清除掉流氓软件、弹窗广告
- 节点偏差Junction Deviation
- 最大信息系数(MIC)
- Flask后端实践 连载十八 Flask输出PDF报表
- XWPFParagraph设置样式
- unity中请问如何点击一下image(image加了button项)变红色再点击一下按钮变成绿色。一直这么循环变色?急求,望大佬指点emmm。c#代码
- [存]MySpace兴衰沉浮启示录:方向混乱 技术欠债多
热门文章
- 升级LTS长期支持版|奇点云数据云平台发布DataSimba R3.8
- Github(Life Restart作者)新作:生火间 网址
- word或excel打开很慢的处理办法
- 群晖NAS使用Docker安装迅雷离线下载出现the active key is not valid.
- Android Studio 之万恶 gradle
- 原生对象、内置对象、宿主对象的区别
- 初学键盘计算机输入时注意,打字练习说明.doc
- Mapper 与 Reducer 解析
- JavaMail实现邮件的发送
- 从零开始设计RISC-V处理器——五级流水线之控制冒险