LU分解求逆

文章轉自：

https://www.cnblogs.com/bigmonkey/p/9555710.html

https://blog.csdn.net/xx_123_1_rj/article/details/39553809

什么是LU分解

在線性代數中， LU分解是矩陣分解的一種，可以將一個矩陣分解為一個單位下三角矩陣和一個上三角矩陣的乘積（有時是它們和一個置換矩陣的乘積）

如果有一個矩陣A，將A表示成下三角矩陣L和上三角矩陣U的乘積，稱為A的LU分解。

 

　　更進一步，我們希望下三角矩陣的對角元素都為1：

LU分解的步驟

　　上一章講到，對於滿秩矩陣A來說，通過左乘一個消元矩陣，可以得到一個上三角矩陣U。

　　可以看到，L實際上就是消元矩陣的逆。容易知道二階矩陣的逆：

 

　　現在假設A是一個3×3矩陣，在不考慮行交換的情況下，通過消元得到上三角矩陣的過程是：

 

 

LU 分解的前提

　　並非所有矩陣都能進行LU分解，能夠LU分解的矩陣需要滿足以下三個條件：

1. 矩陣是方陣（LU分解主要是針對方陣）；
2. 矩陣是可逆的，也就是該矩陣是滿秩矩陣，每一行都是獨立向量；
3. 消元過程中沒有0主元出現，也就是消元過程中不能出現行交換的初等變換

 示例

　　

　　如果A存在LU分解存，a,b滿足什么條件？

　　使用消元法逐一消去主元：

　　由於E₃₁ 中出現了 –b/a，所以a ≠ 0

 

　　b可以是任意常數。

 

具體的算法流程可以是：

 

1)進行LU分解；

2)對分解后的L陣（下三角矩陣）和U陣（上三角矩陣）進行求逆;

3)L陣的逆矩陣和U陣的逆矩陣相乘，即可求得原來矩陣的逆。即：

程序實現如下：

 

#include<stdio.h>
#include <string.h>
#define N 4
void main()
{ 
    float a[N][N];
    float L[N][N],U[N][N],out[N][N], out1[N][N];
    float r[N][N],u[N][N];
    memset( a , 0 , sizeof(a));  //初始化
    memset( L , 0 , sizeof(L));
    memset( U , 0 , sizeof(U));
    memset( r , 0 , sizeof(r));
    memset( u , 0 , sizeof(u));
    int n=N;
    int k,i,j;
    int flag=1;
    float s,t;
    ////////////////////input a matrix////
    printf("input A=\n");
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            scanf("%f",&a[i][j]);
        //////////////////figure the input matrix//////////////////////////
        printf("輸入矩陣:\n");
        for(i=0;i<n;i++)
        {
            for (j = 0; j < n; j++)
            {
                printf("%lf ", a[i][j]);
            }
            printf("\n");
        }
        for(j=0;j<n;j++)
            a[0][j]=a[0][j];  //計算U矩陣的第一行
        
        for(i=1;i<n;i++)
            a[i][0]=a[i][0]/a[0][0];   //計算L矩陣的第1列
        for(k=1;k<n;k++)
        {
            for(j=k;j<n;j++)
            {
                s=0;
                for (i=0;i<k;i++)
                    s=s+a[k][i]*a[i][j];   //累加
                a[k][j]=a[k][j]-s; //計算U矩陣的其他元素
            }
            for(i=k+1;i<n;i++)
            {
                t=0;
                for(j=0;j<k;j++)
                    t=t+a[i][j]*a[j][k];   //累加
                a[i][k]=(a[i][k]-t)/a[k][k];    //計算L矩陣的其他元素
            }
        }
        for(i=0;i<n;i++)
            for(j=0;j<n;j++)
            { 
                if(i>j)
                { 
                    L[i][j]=a[i][j]; 
                    U[i][j]=0;
                }//如果i>j，說明行大於列，計算矩陣的下三角部分，得出L的值，U的//為0
                else
                {
                    U[i][j]=a[i][j];
                    if(i==j) 
                        L[i][j]=1;  //否則如果i<j，說明行小於列，計算矩陣的上三角部分，得出U的//值，L的為0
                    else 
                        L[i][j]=0;
                }
            }  
            if(U[1][1]*U[2][2]*U[3][3]*U[4][4]==0)
            {
                flag=0;
                printf("\n逆矩陣不存在");}
            if(flag==1)
            {
                /////////////////////求L和U矩陣的逆
                for (i=0;i<n;i++) /*求矩陣U的逆 */
                {
                    u[i][i]=1/U[i][i];//對角元素的值，直接取倒數
                    for (k=i-1;k>=0;k--)
                    {
                        s=0;
                        for (j=k+1;j<=i;j++)
                            s=s+U[k][j]*u[j][i];
                        u[k][i]=-s/U[k][k];//迭代計算，按列倒序依次得到每一個值，
                    }
                }
                for (i=0;i<n;i++) //求矩陣L的逆 
                {
                    r[i][i]=1; //對角元素的值，直接取倒數，這里為1
                    for (k=i+1;k<n;k++)
                    {
                        for (j=i;j<=k-1;j++)
                            r[k][i]=r[k][i]-L[k][j]*r[j][i];   //迭代計算，按列順序依次得到每一個值
                    }
                }
                /////////////////繪制矩陣LU分解后的L和U矩陣///////////////////////
                printf("\nLU分解后L矩陣:");
                for(i=0;i<n;i++)
                { 
                    printf("\n");
                    for(j=0;j<n;j++)
                        printf(" %lf",L[i][j]);
                }
                printf("\nLU分解后U矩陣:");
                for(i=0;i<n;i++)
                { 
                    printf("\n");
                    for(j=0;j<n;j++)
                        printf(" %lf",U[i][j]);
                }
                printf("\n");
                ////////繪制L和U矩陣的逆矩陣
                printf("\nL矩陣的逆矩陣:");
                for(i=0;i<n;i++)
                { 
                    printf("\n");
                    for(j=0;j<n;j++)
                        printf(" %lf",r[i][j]);
                }
                printf("\nU矩陣的逆矩陣:");
                for(i=0;i<n;i++)
                { 
                    printf("\n");
                    for(j=0;j<n;j++)
                        printf(" %lf",u[i][j]);
                }
                printf("\n");
                //驗證將L和U相乘，得到原矩陣
                printf("\nL矩陣和U矩陣乘積\n");
                for(i=0;i<n;i++)
                {
                    for(j=0;j<n;j++)
                    {
                        out[i][j]=0;
                    }
                }
                for(i=0;i<n;i++)
                {
                    for(j=0;j<n;j++)
                    {
                        for(k=0;k<n;k++)
                        {
                            out[i][j]+=L[i][k]*U[k][j];
                        }
                    }  
                }
                for(i=0;i<n;i++)
                {
                    for(j=0;j<n;j++)
                    {
                        printf("%lf\t",out[i][j]);
                    }
                    printf("\r\n");
                }
                //////////將r和u相乘，得到逆矩陣
                printf("\n原矩陣的逆矩陣:\n");
                for(i=0;i<n;i++)
                {
                    for(j=0;j<n;j++)
                    {out1[i][j]=0;}
                }
                for(i=0;i<n;i++)
                {
                    for(j=0;j<n;j++)
                    {
                        for(k=0;k<n;k++)            
                        {
                            out1[i][j]+=u[i][k]*r[k][j];
                        }
                    }  
                }
                for(i=0;i<n;i++)
                {
                    for(j=0;j<n;j++)
                    {
                        printf("%lf\t",out1[i][j]);
                    }
                    printf("\r\n");
                }
            }
 }