矩陣求逆的幾種方法總結（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==1)
    {
        return arcs[0][0];
    }
    double ans = 0;
    double temp[N][N]={0.0};
    int i,j,k;
    for(i=0;i<n;i++)
    {
        for(j=0;j<n-1;j++)
        {
            for(k=0;k<n-1;k++)
            {
                temp[j][k] = arcs[j+1][(k>=i)?k+1:k];
                
            }
        }
        double t = getA(temp,n-1);
        if(i%2==0)
        {
            ans += arcs[0][i]*t;
        }
        else
        {
            ans -=  arcs[0][i]*t;
        }
    }
    return ans;
}

//計算每一行每一列的每個元素所對應的余子式，組成A*
void  getAStart(double arcs[N][N],int n,double ans[N][N])
{
    if(n==1)
    {
        ans[0][0] = 1;
        return;
    }
    int i,j,k,t;
    double temp[N][N];
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
        {
            for(k=0;k<n-1;k++)
            {
                for(t=0;t<n-1;t++)
                {
                    temp[k][t] = arcs[k>=i?k+1:k][t>=j?t+1:t];
                }
            }
            
            
            ans[j][i]  =  getA(temp,n-1);  //此處順便進行了轉置
            if((i+j)%2 == 1)
            {
                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(0==flag)
    {
        cout<< "原矩陣行列式為0，無法求逆。請重新運行" <<endl;
        return false;//如果算出矩陣的行列式為0，則不往下進行
    }
    else
    {
        getAStart(src,n,t);
        for(int i=0;i<n;i++)
        {
            for(int j=0;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(0));
    for(int i=0; i<N ;i++)
    {
        for(int j=0; j<N;j++)
        {
            matrix_before[i][j]=rand()%100 *0.01;
        }
    }
    
    cout<<"原矩陣："<<endl;
    
    for(int i=0; i<N ;i++)
    {
        for(int j=0; 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 0;
    
    
    cout<<"逆矩陣："<<endl;
    
    for(int i=0; i<row ;i++)
    {
        for(int j=0; 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=0; i<N ;i++)
    {
        for(int j=0; j<N;j++)
        {
            //cout << matrix_before[i][j] <<" ";
            cout << *(*(matrix_before+i)+j)<<" ";
        }
        cout<<endl;
    }
    
    
    return 0;
}

 

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=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            for(int k=0;k<N;k++)
            {
                C[i*N+j] += A[i*N+k]*B[k*N+j];
            }
        }
    }

    //若絕對值小於10^-10,則置為0（這是我自己定的）
    for(int i=0;i<N*N;i++)
    {
        if(abs(C[i])<pow(10,-10))
        {
            C[i]=0;
        }
    }

    return C;
}

//LUP分解
void LUP_Descomposition(double A[N*N],double L[N*N],double U[N*N],int P[N])
{
    int row=0;
    for(int i=0;i<N;i++)
    {
        P[i]=i;
    }
    for(int i=0;i<N-1;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(0==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=0;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+1;j<N;j++)
        {
            l=A[j*N+i]/u;
            A[j*N+i]=l;
            for(int k=i+1;k<N;k++)
            {
                A[j*N+k]=A[j*N+k]-A[i*N+k]*l;
            }
        }

    }

    //構造L和U
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<=i;j++)
        {
            if(i!=j)
            {
                L[i*N+j]=A[i*N+j];
            }
            else
            {
                L[i*N+j]=1;
            }
        }
        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 = 0;i < N;i++)
    {
        y[i] = b[P[i]];
        for(int j = 0;j < i;j++)
        {
            y[i] = y[i] - L[i*N+j]*y[j];
        }
    }
    //反向替換
    for(int i = N-1;i >= 0; i--)
    {
        x[i]=y[i];
        for(int j = N-1;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=0; 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=0;i<N;i++)
    {
        double *L=new double[N*N]();
        double *U=new double[N*N]();
        int *P=new int[N]();

        //構造單位陣的每一列
        for(int i=0;i<N;i++)
        {
            b[i]=0;
        }
        b[i]=1;

        //每次都需要重新將A復制一份
        for(int i=0;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(0));
    for(int i=0; i<N ;i++)
    {
        for(int j=0; j<N;j++)
        {
            A[i*N+j]=rand()%100 *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=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            cout<< A[i*N+j]<< " " ;
        }
        cout<<endl;
    }

    cout<< "inv_A:" <<endl;
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            cout<< invOfA[i*N+j]<< " " ;
        }
        cout<<endl;
    }

    cout<< "E_test:" <<endl;    
    for(int i=0;i<N;i++)
    {
        for(int j=0;j<N;j++)
        {
            cout<< E_test[i*N+j]<< " " ;
        }
        cout<<endl;
    }

    return 0;
}

 

 

兩種方法運行時間測試樣例（運行環境不同可能會有不同結果，我的主頻是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