矩阵求逆的几种方法总结（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++程序：

 #include <iostream>
 #include <ctime>    //用于产生随机数据的种子 #define N 3    //测试矩阵维数定义 //按第一行展开计算|A|
 double getA(double arcs[N][N],int n)
 {
     if(n==)
     {
         return arcs[][];
     }
     double ans = ;
     double temp[N][N]={0.0};
     int i,j,k;
     for(i=;i<n;i++)
     {
         for(j=;j<n-;j++)
         {
             for(k=;k<n-;k++)
             {
                 temp[j][k] = arcs[j+][(k>=i)?k+:k];             }
         }
         double t = getA(temp,n-);
         if(i%==)
         {
             ans += arcs[][i]*t;
         }
         else
         {
             ans -=  arcs[][i]*t;
         }
     }
     return ans;
 } //计算每一行每一列的每个元素所对应的余子式，组成A*
 void  getAStart(double arcs[N][N],int n,double ans[N][N])
 {
     if(n==)
     {
         ans[][] = ;
         return;
     }
     int i,j,k,t;
     double temp[N][N];
     for(i=;i<n;i++)
     {
         for(j=;j<n;j++)
         {
             for(k=;k<n-;k++)
             {
                 for(t=;t<n-;t++)
                 {
                     temp[k][t] = arcs[k>=i?k+:k][t>=j?t+:t];
                 }
             }             ans[j][i]  =  getA(temp,n-);  //此处顺便进行了转置
             if((i+j)% == )
             {
                 ans[j][i] = - ans[j][i];
             }
         }
     }
 } //得到给定矩阵src的逆矩阵保存到des中。
 bool GetMatrixInverse(double src[N][N],int n,double des[N][N])
 {
     double flag=getA(src,n);
     double t[N][N];
     if(==flag)
     {
         cout<< "原矩阵行列式为0，无法求逆。请重新运行" <<endl;
         return false;//如果算出矩阵的行列式为0，则不往下进行
     }
     else
     {
         getAStart(src,n,t);
         for(int i=;i<n;i++)
         {
             for(int j=;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[N][N]{};//{1,2,3,4,5,6,7,8,9};     //随机数据，可替换
     srand((unsigned)time());
     for(int i=; i<N ;i++)
     {
         for(int j=; j<N;j++)
         {
             matrix_before[i][j]=rand()% *0.01;
         }
     }     cout<<"原矩阵："<<endl;     for(int i=; i<N ;i++)
     {
         for(int j=; j<N;j++)
         {
             //cout << matrix_before[i][j] <<" ";
             cout << *(*(matrix_before+i)+j)<<" ";
         }
         cout<<endl;
     }     double matrix_after[N][N]{};
     flag=GetMatrixInverse(matrix_before,N,matrix_after);
     if(false==flag)
         return ;     cout<<"逆矩阵："<<endl;     for(int i=; i<row ;i++)
     {
         for(int j=; 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=; i<N ;i++)
     {
         for(int j=; j<N;j++)
         {
             //cout << matrix_before[i][j] <<" ";
             cout << *(*(matrix_before+i)+j)<<" ";
         }
         cout<<endl;
     }     return ;
 }

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=;i<N;i++)
     {
         for(int j=;j<N;j++)
         {
             for(int k=;k<N;k++)
             {
                 C[i*N+j] += A[i*N+k]*B[k*N+j];
             }
         }
     }     //若绝对值小于10^-10,则置为0（这是我自己定的）
     for(int i=;i<N*N;i++)
     {
         if(abs(C[i])<pow(,-))
         {
             C[i]=;
         }
     }     return C;
 } //LUP分解
 void LUP_Descomposition(double A[N*N],double L[N*N],double U[N*N],int P[N])
 {
     int row=;
     for(int i=;i<N;i++)
     {
         P[i]=i;
     }
     for(int i=;i<N-;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(==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=;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+;j<N;j++)
         {
             l=A[j*N+i]/u;
             A[j*N+i]=l;
             for(int k=i+;k<N;k++)
             {
                 A[j*N+k]=A[j*N+k]-A[i*N+k]*l;
             }
         }     }     //构造L和U
     for(int i=;i<N;i++)
     {
         for(int j=;j<=i;j++)
         {
             if(i!=j)
             {
                 L[i*N+j]=A[i*N+j];
             }
             else
             {
                 L[i*N+j]=;
             }
         }
         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[N],double b[N])
 {
     double *x=new double[N]();
     double *y=new double[N]();     //正向替换
     for(int i = ;i < N;i++)
     {
         y[i] = b[P[i]];
         for(int j = ;j < i;j++)
         {
             y[i] = y[i] - L[i*N+j]*y[j];
         }
     }
     //反向替换
     for(int i = N-;i >= ; i--)
     {
         x[i]=y[i];
         for(int j = N-;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=; 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[N]();//矩阵逆的各列
     //double *B    =new double[N*N]();
     double *b    =new double[N]();//b阵为B阵的列矩阵分量     for(int i=;i<N;i++)
     {
         double *L=new double[N*N]();
         double *U=new double[N*N]();
         int *P=new int[N]();         //构造单位阵的每一列
         for(int i=;i<N;i++)
         {
             b[i]=;
         }
         b[i]=;         //每次都需要重新将A复制一份
         for(int i=;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());
     for(int i=; i<N ;i++)
     {
         for(int j=; j<N;j++)
         {
             A[i*N+j]=rand()% *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=;i<N;i++)
     {
         for(int j=;j<N;j++)
         {
             cout<< A[i*N+j]<< " " ;
         }
         cout<<endl;
     }     cout<< "inv_A:" <<endl;
     for(int i=;i<N;i++)
     {
         for(int j=;j<N;j++)
         {
             cout<< invOfA[i*N+j]<< " " ;
         }
         cout<<endl;
     }     cout<< "E_test:" <<endl;
     for(int i=;i<N;i++)
     {
         for(int j=;j<N;j++)
         {
             cout<< E_test[i*N+j]<< " " ;
         }
         cout<<endl;
     }     return ;
 }

两种方法运行时间测试样例（运行环境不同可能会有不同结果，我的主频是2.6GHz,内存8g。时间单位：毫秒ms）

个人认为LU分解法的两个1ms其实是不准确的（实际应远小于1ms，有兴趣可以试试看）。

三种方法的复杂度分析：

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