求逆矩阵（C++）

求矩阵的逆常见的一般有三种方法（考研常见）：待定系数法、高斯-约旦消元法和伴随矩阵求逆矩阵。

1. 待定系数法：
   假设矩阵A：
   $$\left[\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\\ 7& 8& 9\end{array}\right]$$
   设其逆矩阵为$A^{-1}$：
   $$\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$$
   则有如下方程组：
   $$\{\begin{array}{l}1a_{11}+2a_{21}+3a_{31}=1\\ 1a_{12}+2a_{22}+3a_{32}=0\\ 1a_{13}+2a_{23}+3a_{33}=0\\ 4a_{11}+5a_{21}+6a_{31}=0\\ 4a_{12}+5a_{22}+6a_{32}=1\\ 4a_{13}+5a_{23}+6a_{33}=0\\ 7a_{11}+8a_{21}+9a_{31}=0\\ 7a_{12}+8a_{22}+9a_{32}=0\\ 7a_{13}+8a_{23}+9a_{33}=1\end{array}$$
   解出该方程组即可得到逆矩阵$A^{-1}$。

2. 伴随矩阵求逆矩阵：
   设矩阵$A$的伴随矩阵为$A^{\ast }$，则有
   $$A^{-1}=\frac{1}{|A|}{A}^{\ast }$$
   其中$|A|$为矩阵$A$的行列式。
   伴随矩阵的概念如下（注意里面元素的下标）：
   $$A^{\ast }=\left[\begin{array}{cccc}{A}_{11}& A_{21}& \cdots & A_{n1}\\ A_{12}& A_{22}& \cdots & A_{n2}\\ ⋮& ⋮& & ⋮\\ A_{1n}& A_{2n}& \cdots & A_{nn}\end{array}\right]$$
   其中$A_{ij}$为矩阵$A$的第i,j行代数余子式。$M_{ij}$余子式为矩阵$A$去除第i行和j列元素后的矩阵，$A_{ij}=(-1{)}^{i+j}{M}_{ij}$。
   因此如果需要使用这个方法求逆矩阵，则先要实现行列式求值，使用伴随矩阵求逆矩阵的代码如下：

   //使用伴随矩阵法求矩阵的逆
   #include <iostream>
   //打印矩阵
   void print_matrix(double **A,int m,int n)
   {
       for (int i = 0; i < m;i++)
       {
           for (int j = 0; j < n;j++)
           {
               std::cout<<A[i][j]<<" ";
           }
           std::cout << std::endl;
       }
   }
   //使用代数余子式求矩阵的值
   double determinant_value(double **D,int n)
   {
       //递归终点
       if(n==1)
       {
           return  D[0][0];
       }
       else if(n==2)
       {
           return D[1][1]*D[0][0]-D[0][1]*D[1][0];
       }
       else{
           double D_value=0;
           for(int k=0;k<n;k++)
           {
               double **D_temp=new double *[n-1];
               int i2=1;//由于是根据第0行第j列展开，所以原本的行列式直接从第一行开始拷贝
               for(int i1=0;i1<n-1;i1++)
               {
                   D_temp[i1]=new double[n-1];
                   int j2=0;
                   for(int j1=0;j1<n-1;j1++)
                   {
                       if(j2==k)
                       {
                           j2++;
                       }//去除第j列
                       D_temp[i1][j1]=D[i2][j2++];
                   }
                   i2++;
               }
               D_value+=((k&0x0001)?(-1):1)*D[0][k]*determinant_value(D_temp,n-1);//计算代数余子式与对应项的乘积
           }
           return D_value;
       }
   }
   //求伴随矩阵 
   void Adjoint_matrices(double **A,double **A_adj,int n)
   {
       double **A_AlgebraicRemainder=new double *[n-1];//暂时存放代数余子式的矩阵
       for (int i = 0; i < n-1;i++)
       {
           A_AlgebraicRemainder[i]=new double[n-1];//初始化分配空间
       }
       for (int x = 0; x < n; x++)
       {
           for (int y = 0; y < n; y++)
           {
               int i1 = 0;
               for (int i = 0; i < n-1; i++)
               {
                   if(i1==x)
                   {
                       i1++;
                   }
                   int j1 = 0;
                   for (int j = 0; j < n-1; j++)
                   {
                       if(j1==y)
                       {
                           j1++;
                       }
                       A_AlgebraicRemainder[i][j]=A[i1][j1++];
                   }
                   i1++;
               }
               //实现了转置
               A_adj[y][x]=((x+y)&0x0001?(-1):1)*determinant_value(A_AlgebraicRemainder, n - 1);
           }
       }
       for(int i=0;i<n-1;i++)
       {
           delete[] A_AlgebraicRemainder[i];
       }
       delete[] A_AlgebraicRemainder;
   }
   
   //求逆矩阵
   void matrix_inverse(double **A,double **A_inverse,int n)
   {
       double A_D=determinant_value(A,n);
       
       if (A_D == 0.0f)
           {
               std::cout << "该矩阵不可逆" << std::endl;
               return;
           }
           else{
               if(n==1)
               {
                   A_inverse[0][0]=1/A[0][0];
                   return;
               }
               
               double **A_adj=new double *[n];
               for (int i = 0; i < n;i++)
               {
                   A_adj[i]=new double[n];
               }
               Adjoint_matrices(A, A_adj, n);
               for (int i = 0; i < n;i++)
               {
                   for (int j = 0; j < n;j++)
                   {
                       A_inverse[i][j]=A_adj[i][j]/A_D;
                   }
               }
               for (int i = 0; i < n; i++)
               {
                   delete[] A_adj[i];
               }
               delete[] A_adj;
       }
   }
   
   int main() 
   {
       int n;
       std::cout<<"Enter the matrix dimension A: ";
       std::cin>>n;//输入数组维度
       double **A=new double *[n];
       std::cout<<"Enter the coefficient matrix:"<<std::endl;
       for(int i=0;i<n;i++)
       {
           A[i]=new double[n];
           for(int j=0;j<n;j++)
           {
               std::cin>>A[i][j];//每次输入一个数字都用空格隔开,输入样例
               //1 2 3\enter
               //4 5 6\enter
               //7 8 9\enter
           }
       }
       double **A_inverse=new double *[n];
       for(int i=0;i<n;i++)
       {
           A_inverse[i]=new double[n];
       }
       matrix_inverse(A, A_inverse, n);
       std::cout<<"The inverse matrix is:"<<std::endl;
       print_matrix(A_inverse, n, n);
       return 0;
   }
   

   

3. 高斯消元法求逆矩阵
   这个方法用于手算比较多，但是代码实现也比较简单，其原理如下：
   设矩阵$A$:
   $$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$$
   在其右侧增加一个单位矩阵$E$，其实成为增广矩阵$\hat{A}=[A|E]$:
   $$\hat{A}=\left[\begin{array}{ccccccc}{a}_{11}& a_{12}& a_{13}& |& 1& 0& 0\\ a_{21}& a_{22}& a_{23}& |& 0& 1& 0\\ a_{31}& a_{32}& a_{33}& |& 0& 0& 1\end{array}\right]$$
   然后通过初等行变换使得增广矩阵左边变为单位矩阵，右侧就变为了$A$的逆矩阵。
   对应的代码如下：

   //这里展示如何用C++求矩阵的逆矩阵
   #include <iostream>
   //矩阵交换某两行
   void matrix_swap_row(double **A,int i,int j,int n)
   {
       double temp;
       for(int k=0;k<n;k++)
       {
           temp=A[i][k];
           A[i][k]=A[j][k];
           A[j][k]=temp;
       }
   }
   //矩阵第i行=矩阵第i行-矩阵第j行*a
   void matrix_minus_inner(double **A,double a,int i,int j,int n)
   {
       for(int k=0;k<n;k++)
       {
           A[i][k]-=a*A[j][k];
       }
   }
   //矩阵求逆
   void matrix_inverse(double **A,double **A_inverse,int n)
   {
       double **A_E=new double*[2*n];
       //构建增广矩阵
       for(int i=0;i<n;i++)
       {
           A_E[i]=new double [n*2];
           for(int j=0;j<n*2;j++)
           {
               if(j<n)
               {
                   A_E[i][j]=A[i][j];
               }
               else if((j-n)==i){
                   A_E[i][j]=1;
               }
               else{
                   A_E[i][j]=0;
               }
           }
       }
       //首先将矩阵化为上三角矩阵
       for(int i=0;i<n;i++)
       {
           if(A_E[i][i]==0)
           {
               for(int k=i+1;k<n;k++)
               {
                   if(A_E[k][i]!=0)
                   {
                       matrix_swap_row(A_E,i,k,n*2);
                       break;
                   }
               }
           }
           for(int j=i+1;j<n;j++)
           {
               matrix_minus_inner(A_E,A_E[j][i]/A_E[i][i],j,i,2*n);
           }
       }
       //打印上三角矩阵进行检验
       // for(int i=0;i<n;i++)
       // {
       //     for(int j=0;j<n*2;j++)
       //     {
       //         std::cout<<A_E[i][j]<<" ";
       //     }
       //     std::cout<<std::endl;
       // }
       //判断矩阵是否可逆
       for(int i=0;i<n;i++)
       {
           if(A_E[i][i]==0)
           {
               std::cout<<"矩阵不可逆"<<std::endl;
               exit(0);
           }
       }
       //将上三角转换为对角矩阵
       for(int j=1;j<n;j++)
       {
           for(int i=0;i<j;i++)
           {
               matrix_minus_inner(A_E,A_E[i][j]/A_E[j][j],i,j,2*n);
           }
       }
       for(int i=0;i<n;i++)
       {
           for(int j=n;j<2*n;j++)
           {
               A_inverse[i][j-n]=A_E[i][j]/A_E[i][i];
           }
       }
       // for(int i=0;i<n;i++)
       // {
       //     for(int j=0;j<n;j++)
       //     {
       //         std::cout<<A_inverse[i][j]<<" ";
       //     }
       //     std::cout<<std::endl;
       // }
   
   }
   int main()
   {
       int n;
       std::cout<<"Enter the matrix dimension A: ";
       std::cin>>n;//输入数组维度
       double **A=new double *[n];
       std::cout<<"Enter the coefficient matrix:"<<std::endl;
       for(int i=0;i<n;i++)
       {
           A[i]=new double[n];
           for(int j=0;j<n;j++)
           {
               std::cin>>A[i][j];//每次输入一个数字都用空格隔开,输入样例
               //1 2 3\enter
               //4 5 6\enter
               //7 8 9\enter
           }
       }
       double **A_inverse=new double *[n];
       for(int i=0;i<n;i++)
       {
           A_inverse[i]=new double[n];
       }
       matrix_inverse(A, A_inverse, n);
       std::cout<<"The inverse matrix of A is:"<<std::endl;
       for(int i=0;i<n;i++)
       {
           for(int j=0;j<n;j++)
           {
               std::cout<<A_inverse[i][j]<<" ";
           }
           std::cout<<std::endl;
       }
       return 0;
   }
   

   运行结果如下：
   

   上面两个结果都与matlab的输出结果一致，说明程序正确。