求逆矩阵 —— LU分解法

LU分解：

算法步骤：

1. 将 A 矩阵分解为 L 下三角矩阵和 U 上三角矩阵。

step1. L对角线填充为1

step2. 

step3. 

step4. U是按行迭代计算，L是按列迭代计算，UL交错计算，且U先L一步

    for  k = 1  to  m-1: {

   

}

2. 分别对 L 和 U 求逆，得到 Linv 和 Uinv.

step1. 列顺序行顺序, for  j = 0  to  m-1, i = j  to  m-1

step2. 列顺序行倒序, for j = 0 to m-1, i = j to 0

3. Ainv = Linv * Uinv. 

实现代码(java)：

public class TestMain {

	public static void main(String[] args) {
		double[][] A = { { 4, 2, 1, 5 }, { 8, 7, 2, 10 }, { 4, 8, 3, 6 }, { 6, 8, 4, 9 } };
		double[][] U = new double[4][4];
		double[][] L = new double[4][4];
		double[][] Uinv = new double[4][4];
		double[][] Linv = new double[4][4];
		print(A);
		int size = 4;
		// fill 1 at L's diagonal
		for (int i = 0; i < size; i++) {
			L[i][i] = 1.0;
		}
		// calculate the U's first row
		for (int j = 0; j < size; j++) {
			U[0][j] = A[0][j];
		}
		// calculate the L's first column
		for (int i = 1; i < size; i++) {
			L[i][0] = A[i][0] / U[0][0];
		}
		// iterative calculation of other rows and columns
		for (int k = 1; k < size; k++) {
			// the k_th row of U
			for (int j = k; j < size; j++) {
				// sum(L_kt*U_tj),t
				double s = 0.0;
				for (int t = 0; t < k; t++) {
					s += L[k][t] * U[t][j];
				}
				// U_kj = A_kj - s
				U[k][j] = A[k][j] - s;
			}
			// the k_th column of L
			for (int i = k; i < size; i++) {
				// sum(L_it*U_tk),t
				double s = 0.0;
				for (int t = 0; t < k; t++) {
					s += L[i][t] * U[t][k];
				}
				// L_ik = (A_ik - s) / U_kk
				L[i][k] = (A[i][k] - s) / U[k][k];
			}
		}
		
		// calculate L_inv
		for (int j = 0; j < size; j++) {
			for (int i = j; i < size; i++) {
				if (i == j) Linv[i][j] = 1 / L[i][j];
				else if (i < j) Linv[i][j] = 0;
				else {
					double s = 0.0;
					for (int k = j; k < i; k++) {
						s += L[i][k] * Linv[k][j];
					}
					Linv[i][j] = -Linv[j][j] * s;
				}
			}
		}
		
		// calculate U_inv
		for (int j = 0; j < size; j++) {
			for (int i = j; i >= 0; i--) {
				if (i == j) Uinv[i][j] = 1 / U[i][j]; 
				else if (i > j) Uinv[i][j] = 0;
				else {
					double s = 0.0;
					for (int k = i + 1; k <= j; k++) {
						s += U[i][k] * Uinv[k][j];
					}
					Uinv[i][j] = -1 / U[i][i] * s;
				}
			}
		}
		
		double[][] inv = new double[4][4];
		for (int i = 0; i < size; i++) {
			for (int j = 0; j < size; j++) {
				for (int k = 0; k < size; k++) {
					inv[i][j] += Uinv[i][k] * Linv[k][j];
				}
			}
		}
		
		print(L);
		print(U);
		print(Linv);
		print(Uinv);
		print(inv);
	}

	public static void print(double[][] M) {
		for (int i = 0; i < M.length; i++) {
			for (int j = 0; j < M[0].length; j++) {
				System.out.printf("%8.2f ", M[i][j]);
			}
			System.out.println();
		}
		System.out.println();
	}

}

程序输出：

    4.00     2.00     1.00     5.00

    8.00     7.00     2.00    10.00

    4.00     8.00     3.00     6.00

    6.00     8.00     4.00     9.00

 

    1.00     0.00     0.00     0.00

    2.00     1.00     0.00     0.00

    1.00     2.00     1.00     0.00

    1.50     1.67     1.25     1.00

 

    4.00     2.00     1.00     5.00

    0.00     3.00     0.00     0.00

    0.00     0.00     2.00     1.00

    0.00     0.00     0.00     0.25

 

    1.00     0.00     0.00     0.00

   -2.00     1.00     0.00     0.00

    3.00    -2.00     1.00     0.00

   -1.92     0.83    -1.25     1.00

 

    0.25    -0.17    -0.13    -4.50

    0.00     0.33    -0.00    -0.00

    0.00     0.00     0.50    -2.00

    0.00     0.00     0.00     4.00

 

    8.83    -3.67     5.50    -4.50

   -0.67     0.33     0.00     0.00

    5.33    -2.67     3.00    -2.00

   -7.67     3.33    -5.00     4.00