矩阵求逆的几种方法总结(C++)

矩阵求逆运算有多种算法:

  1. 伴随矩阵的思想,分别算出其伴随矩阵和行列式,再算出逆矩阵;
  2. 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,有兴趣可以试试看)。

三种方法的复杂度分析:

  1. 伴随矩阵法:此法的时间复杂度主要来源于计算行列式,由于计算行列式的函数为递归形式,其复杂度为O(n2)[参见这里],而整体算法需要计算每个元素的代数余子式,时间复杂度直接扩大n2倍,变为O(n4)。而递归算法本身是需要占用栈空间的,因此需要注意:当矩阵的维数较大时,随着递归深度的加大,临时占用的空间将会越来越多,甚至可能会出现栈不够用的情况(当然本次实现没有遇到,因为此时的时间开销实在令人难以忍受)!
  2. LU分解法:此法主要是分解过程耗时,求解三角矩阵的时间复杂度是O(n2),分解过程是O(n3),总体来说和高斯消元法差不多,但是避免了高斯消元法的主元素为0的过程。 为了节省空间,A=LU分解的元素存放在A的矩阵中(因为当用过了a[i][j]元素后,便不再用了,所以可以占用原矩阵A的空间)。但是有利就有弊,考虑如果是上千个元素的矩阵,引用传参,这样就改变原矩阵了,因此程序中使用A_mintor作为副本进行使用。另外,可以看出,当矩阵维数超过某值时,内存空间便不够用了(具体是多少没有试验)。还需注意的一点是:程序中未对矩阵是否奇异进行检查,如果矩阵奇异,就不应再进行下去了。
  3. LU分解法中,还可以先分别求出U和L的逆,再相乘,此法其实与常规LU分解法差不多。

其他:

文章中用到了矩阵的原地转置算法,具体请参考第4篇文献,这种方法降低了空间复杂度。

 

需要注意的问题:

本文介绍的方法new了一些指针,未释放,会出现内存泄漏,使用前请释放掉。

 

本文参考了以下几篇文章:

 

posted @ 2017-02-20 20:00 xiaoxi666 阅读(...) 评论(...) 编辑 收藏
TOP