求逆矩阵 —— LU分解法「建议收藏」

全栈程序员-站长 • 2022年8月21日下午5:00 • 未分类 • 阅读 12

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元售后保障童叟无欺

LU分解：

$A=L*U\Rightarrow A^{-1} = (L*U)^{-1}=L^{-1}*U^{-1}$

算法步骤：

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

step1. L对角线填充为1

step2. $U_{0,j}=A_{0,j}(j=0,...,m-1)$

step3. $L_{i,0}=\frac{A_{i,0}}{U_{0,0}},(i=1,...,m-1)$

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

for k = 1 to m-1: {

$U_{k,j}=\frac{A_{k,j}-\sum_{t=0}^{k-1}(L_{i,t}U_{t,j})}{L_{k,k}}=A_{k,j}-\sum_{t=0}^{k-1}(L_{k,t} U_{t,j}),(j=k,...,m-1)$

$L_{i,k}=\frac{A_{i,k}-\sum_{t=0}^{k-1}(L_{i,t}U_{t,k})}{U_{k,k}},(i=k,...,m-1)$

}

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

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

求逆矩阵 —— LU分解法「建议收藏」 step2. 列顺序行倒序, for j = 0 to m-1, i = j to 0

求逆矩阵 —— LU分解法「建议收藏」

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