【图像处理】矩阵运算代码实现2-矩阵求逆

这篇总结是《矩阵运算代码实现1》的后续，主要对矩阵求逆的算法及代码实现进行描述。具体如下。

矩阵求逆算法

这里采用的是LU分解对矩阵进行求逆。原理如下：

inv(A)=inv(LU)= inv(U)inv(L)
将原矩阵A分解成两个三角矩阵，上三角矩阵L和下三角矩阵U。通过求解U和L对应的逆矩阵，即可求得相应A的逆矩阵。

算法分成以下几个步骤：
1）矩阵的LU分解

对于LU上每个位置的值，可以用以下公式进行计算得到：

2）求对应LU的逆矩阵

需要注意的是，在求uji的时候，是采用逆序的方式，如下图所示，

求解过程：【先求u12，再求u02】，具体可参照实现的代码。

3）求A的逆矩阵

在求得L 和 U 对应的逆矩阵之后，可以得到对应的A的逆矩阵

该算法主要参考自：http://www.doc88.com/p-8498050172649.html
但该文中存在部分错漏，这里已经修改，最终描述如上。

具体实现，可参考如下代码：

代码实现

int matInv(double *a,int ra,double *inv){
    if(a==NULL||inv==NULL){
        return 0;
    }
    int i,j,k,r;
    double *L = (double*)malloc(sizeof(double)*ra*ra);if(L==NULL)return 0;
    double *U = (double*)malloc(sizeof(double)*ra*ra);if(U==NULL)return 0;
    double *l = (double*)malloc(sizeof(double)*ra*ra);if(l==NULL)return 0;
    double *u = (double*)malloc(sizeof(double)*ra*ra);if(u==NULL)return 0;

    double sum;

    // L U 初始化
    for(i=0;i<ra*ra;i++){
        *(L+i)=0.0;
        *(U+i)=0.0;
        *(l+i)=0.0;
        *(u+i)=0.0;
    }
    // L的对角 U的第一行
    for(i=0;i<ra;i++){
        *(L+i*ra+i)=1;
        *(U+i)=*(a+i);
    }

    // L的第一列
    for(i=1;i<ra;i++){
        *(L+ra*i)=*(a+ra*i)/(*U);
    }

    // U的第二行 L的第二列，... 
    for(r=1;r<ra;r++){
        for(j=r;j<ra;j++){
            sum=0.0;
            for(k=0;k<r;k++){
                sum+=*(L+r*ra+k)*(*(U+k*ra+j));
            }
            *(U+r*ra+j)=*(a+ra*r+j)-sum;
        }
        for(i=r+1;i<ra;i++){
            sum=0.0;
            for(k=0;k<r;k++){
                sum+=*(L+i*ra+k)*(*(U+k*ra+r));
            }
            *(L+i*ra+r)=(*(a+i*ra+r)-sum)/(*(U+r*ra+r));
        }
    }

//  printf("L:\r\n");
//  matPrint(L,ra,ra);
//  printf("U:\r\n");
//  matPrint(U,ra,ra);

    // 求U L 的逆 u l
    for(i=0;i<ra;i++){
        *(l+i*ra+i)=1/(*(L+i*ra+i));
        *(u+i*ra+i)=1/(*(U+i*ra+i));
    }
    for(j=0;j<ra;j++){
        for(i=0;i<j;i++){
            sum = 0.0;
            for(k=i;k<j;k++){
                sum+=*(L+j*ra+k)*(*(l+k*ra+i));
            }
            *(l+j*ra+i)=-(*(l+i*ra+i))*sum;
        }       
    }
//  printf("l:\r\n");
//  matPrint(l,ra,ra);
    for(i=1;i<ra;i++){
        for(j=i-1;j>=0;j--){
            sum = 0.0;
            for(k=j+1;k<=i;k++){
                sum+=*(U+j*ra+k)*(*(u+k*ra+i));
            }
            *(u+j*ra+i)=-sum*(*(u+j*ra+j));
        }
    }

    printf("u:\r\n");
    matPrint(u,ra,ra);

    matMul(u,l,ra,ra,ra,ra,inv);

    free(L);
    free(U);
    free(l);
    free(u);


}