矩阵求逆的几种方法总结（C++）
矩阵求逆的几种方法总结（C++）
xiaoxi666
发布于 2018-10-29 17:11:14
发布于 2018-10-29 17:11:14
矩阵求逆运算有多种算法：
- 伴随矩阵的思想，分别算出其伴随矩阵和行列式，再算出逆矩阵；
- 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
本文参与 腾讯云自媒体分享计划，分享自作者个人站点/博客。
原始发表：2017-02-20 ，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读