LU分解的矩阵逆运算

算法名称：矩阵求逆（基于LU分解法）

 

LU分解算法评价：

       LU分解大约需要执行N³/3次内层循环（每次包括一次乘法和一次加法）。这是求解一个（或少量几个）右端项时的运算次数，它要比Gauss-Jordan消去法快三倍，比不计算逆矩阵的Gauss-Jordan法快1.5倍。

       当要求解逆矩阵时，总的运算次数（包括向前替代和回代部分）为N³，与Gauss-Jordan法相同。

 

算法描述：

简言之，我们只需对原始矩阵进行一次LU分解，然后变换右端向量b就可以了，即设我们的原始矩阵为4×4阶方阵，那么我们的b依次取

然后重新排列成的矩阵就是逆矩阵了。

 

运行示例：

Origin matrix:

 | 0.0 2.0 0.0 1.0 |

 | 2.0 2.0 3.0 2.0 |

 | 4.0 -3.0 0.0 1.0 |

 | 6.0 1.0 -6.0 -5.0 |

-----------------------------------------------

Its inverse matrix:

 | -0.025641025641025623 0.1282051282051282 0.08974358974358977 0.0641025641025641 |

 | 0.17948717948717946 0.10256410256410259 -0.12820512820512822 0.05128205128205129 |

 | -0.5299145299145299 0.3162393162393163 -0.14529914529914528 -0.00854700854700854 |

 | 0.6410256410256411 -0.20512820512820518 0.25641025641025644 -0.10256410256410257 |

-----------------------------------------------

示例程序：

package  com.nc4nr.chapter02.matrixinver;

public   class  MatrixInver  {

    double[][] a = {
            {0.0, 2.0, 0.0, 1.0},
            {2.0, 2.0, 3.0, 2.0},
            {4.0, -3.0, 0.0, 1.0},
            {6.0, 1.0, -6.0, -5.0}
    };
    
    double[] b = null;
    
    int anrow = 4;
    double[] vv = new double[anrow];
    int[] indx = new int[anrow];
    
    private void lucmp() {
        int n = anrow, imax = 0;
        for (int i = 0; i < n; i++) {
            double big = 0.0;
            for (int j = 0; j < n; j++) {
                double temp = Math.abs(a[i][j]);
                if (temp > big) big = temp;
            }
            vv[i] = 1.0 / big;
        }
        for (int j = 0; j < n; j++) {
            for (int i = 0; i < j; i++) {
                double sum = a[i][j];
                for (int k = 0; k < i; k++) sum -= a[i][k] * a[k][j];
                a[i][j] = sum;
            }
            double big = 0.0;
            for (int i = j; i < n; i++) {
                double sum = a[i][j];
                for (int k = 0; k < j; k++) sum -= a[i][k] * a[k][j];
                a[i][j] = sum;
                double dum = vv[i] * Math.abs(sum);
                if (dum >= big) {
                    big = dum;
                    imax = i;
                }
            }
            if (j != imax) {
                for (int i = 0; i < n; i++) {
                    double mid = a[imax][i];
                    a[imax][i] = a[j][i];
                    a[j][i] = mid;
                }
                double mid = vv[j];
                vv[j] = vv[imax];
                vv[imax] = mid;
            }
            indx[j] = imax;
            if (j != n - 1) {
                double dum = 1.0/a[j][j];
                for (int i = j + 1; i < n; i++) a[i][j] *= dum; 
            }
        }
    }
    
    private void lubksb(double[] b) {
        int n = anrow, ii = 0;
        // y
        for (int i = 0; i < n; i++) {    
            int ip = indx[i];        
            double sum = b[ip];
            b[ip] = b[i];
            if (ii != 0)
                for (int j = ii - 1; j < i; j++) sum -= a[i][j] * b[j];
            else
                ii = i + 1;
            b[i] = sum;
        }
        // x
        for (int i = n - 1; i >= 0; i--) {
            double sum = b[i];
            for (int j = i + 1; j < n; j++) sum -= a[i][j]*b[j];
            b[i] = sum / a[i][i];
        }
    }
    
    private void output(double a[][], int anrow) {
        for (int i = 0; i < anrow; i++) {
            System.out.println(" | " + a[i][0] + " " + 
                    a[i][1] + " " + 
                    a[i][2] + " " + 
                    a[i][3] + " | ");
        }
        System.out.println("-----------------------------------------------");
    }
    
    public MatrixInver() {
        System.out.println("Origin matrix:");
        output(a,4);
        lucmp();
        double[] b = new double[anrow];
        double[][] y = new double[anrow][anrow];
        for (int i = 0; i < anrow; i++) { 
            for (int j = 0; j < anrow; j++) b[j] = 0;
            b[i] = 1.0;
            lubksb(b);
            for (int j = 0; j < anrow; j++) y[j][i] = b[j];
        }
        System.out.println("Its inverse matrix:");
        output(y,4);
    }
    
    public static void main(String[] args) {
        new MatrixInver();
    }

}

http://blog.csdn.net/BoyMgl/archive/2007/12/03/1914288.aspx