算法概述
LU分解（LU Decomposition / LU Factorization）是数值线性代数中最基础的矩阵分解技术，由Turing在1948年正式提出。其核心思想是将一个方阵 (A) 分解为一个下三角矩阵 (L)（Lower triangular，对角线元素为1）和一个上三角矩阵 (U)（Upper triangular）的乘积：
[A = LU]
一旦完成LU分解，求解线性方程组 (Ax = b) 就可以转化为两步简单的三角方程组求解：先解 (Ly = b)（前向替换），再解 (Ux = y)（后向替换）。相比于直接高斯消元的 (O(n^3))，LU分解将求解多个右端项的问题优化为一次 (O(n^3)) 分解加多次 (O(n^2)) 求解。
本文用Java实现带部分主元选取（Partial Pivoting）的LU分解，并展示其在方程组求解、矩阵求逆和行列式计算中的应用。
核心原理
基本LU分解（Doolittle分解）
对于 (n \times n) 矩阵 (A)，通过高斯消元将其转化为上三角矩阵 (U)，同时记录消元乘数构成下三角矩阵 (L)。
第 (k) 步消元：
1. 对第 (k) 列，计算乘数 (m_{ik} = a_{ik}^{(k)} / a_{kk}^{(k)})（(i = k+1, …, n)）
2. 对行 (i) 执行：(R_i \leftarrow R_i – m_{ik} \cdot R_k)
3. 将乘数 (m_{ik}) 存入 (L) 的对应位置
部分主元选取（Partial Pivoting）
为避免除数过小导致数值不稳定，每一步在第 (k) 列下方选取绝对值最大的元素作为主元，并交换行。此时分解变为：
[PA = LU]
其中 (P) 是置换矩阵，记录行交换历史。
Java实现
项目结构
lu-decomposition/
├── LUDecomposition.java // LU分解核心实现
├── LinearSystemSolver.java // 线性方程组求解器
├── MatrixOps.java // 矩阵运算工具
└── Main.java // 测试入口
核心实现：LUDecomposition.java
/**
* LU分解实现（带部分主元选取）
* 支持：分解、求解线性方程组、矩阵求逆、行列式计算
*/
public class LUDecomposition {
private final double[][] L; // 下三角矩阵
private final double[][] U; // 上三角矩阵
private final int[] pivot; // 行置换记录
private final int n;
private boolean decomposed;
private double detSign; // 行列式符号（由行交换次数决定）
public LUDecomposition(double[][] A) {
if (A == null || A.length == 0 || A.length != A[0].length) {
throw new IllegalArgumentException("输入必须是方阵");
}
this.n = A.length;
this.L = new double[n][n];
this.U = new double[n][n];
this.pivot = new int[n];
this.decomposed = false;
// 复制A到U，准备分解
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
U[i][j] = A[i][j];
}
pivot[i] = i;
L[i][i] = 1.0; // L的对角线为1
}
}
/**
* 执行Doolittle LU分解（带部分主元）
* 时间复杂度: O(n^3)
*/
public void decompose() {
if (decomposed) return;
detSign = 1.0;
for (int k = 0; k < n; k++) {
// 步骤1：在第k列下方找主元（绝对值最大）
int maxRow = k;
double maxVal = Math.abs(U[k][k]);
for (int i = k + 1; i < n; i++) {
if (Math.abs(U[i][k]) > maxVal) {
maxVal = Math.abs(U[i][k]);
maxRow = i;
}
}
// 步骤2：如果主元不在当前行，交换行
if (maxRow != k) {
swapRows(U, k, maxRow);
// 只交换L中已计算的部分（前k列）
for (int j = 0; j < k; j++) {
double temp = L[k][j];
L[k][j] = L[maxRow][j];
L[maxRow][j] = temp;
}
// 交换置换记录
int tempPivot = pivot[k];
pivot[k] = pivot[maxRow];
pivot[maxRow] = tempPivot;
detSign *= -1; // 行交换改变行列式符号
}
// 检查奇异矩阵
if (Math.abs(U[k][k]) < 1e-12) {
throw new ArithmeticException("矩阵奇异或接近奇异，无法分解");
}
// 步骤3：计算消元乘数并存入L，执行消元更新U
for (int i = k + 1; i < n; i++) {
double multiplier = U[i][k] / U[k][k];
L[i][k] = multiplier;
// 消元：R_i = R_i - multiplier * R_k
for (int j = k; j < n; j++) {
U[i][j] -= multiplier * U[k][j];
}
}
}
decomposed = true;
}
/**
* 求解线性方程组 Ax = b
* 分解后每解一个右端项仅需 O(n^2)
*/
public double[] solve(double[] b) {
if (!decomposed) decompose();
if (b.length != n) {
throw new IllegalArgumentException("右端项维度不匹配");
}
// 步骤1：应用行置换 Pb
double[] pb = new double[n];
for (int i = 0; i < n; i++) {
pb[i] = b[pivot[i]];
}
// 步骤2：前向替换解 Ly = Pb
double[] y = forwardSubstitute(L, pb);
// 步骤3：后向替换解 Ux = y
return backwardSubstitute(U, y);
}
/**
* 求解多右端项 AX = B
*/
public double[][] solveMultiple(double[][] B) {
if (!decomposed) decompose();
if (B.length != n) {
throw new IllegalArgumentException("矩阵维度不匹配");
}
int m = B[0].length;
double[][] X = new double[n][m];
for (int col = 0; col < m; col++) {
double[] b = new double[n];
for (int i = 0; i < n; i++) {
b[i] = B[i][col];
}
double[] x = solve(b);
for (int i = 0; i < n; i++) {
X[i][col] = x[i];
}
}
return X;
}
/**
* 计算矩阵行列式
*/
public double determinant() {
if (!decomposed) decompose();
double det = detSign;
for (int i = 0; i < n; i++) {
det *= U[i][i];
}
return det;
}
/**
* 计算矩阵逆 A^(-1)
* 通过求解 n 个单位向量右端项实现
*/
public double[][] inverse() {
if (!decomposed) decompose();
double[][] inv = new double[n][n];
// 对每一列求解 A * x = e_i
for (int col = 0; col < n; col++) {
double[] e = new double[n];
e[col] = 1.0;
double[] x = solve(e);
for (int i = 0; i < n; i++) {
inv[i][col] = x[i];
}
}
return inv;
}
/**
* 前向替换：解下三角方程组 L*y = b
*/
private double[] forwardSubstitute(double[][] L, double[] b) {
double[] y = new double[n];
for (int i = 0; i < n; i++) {
double sum = b[i];
for (int j = 0; j < i; j++) {
sum -= L[i][j] * y[j];
}
y[i] = sum / L[i][i]; // L[i][i] = 1
}
return y;
}
/**
* 后向替换：解上三角方程组 U*x = y
*/
private double[] backwardSubstitute(double[][] U, double[] y) {
double[] x = new double[n];
for (int i = n - 1; i >= 0; i--) {
double sum = y[i];
for (int j = i + 1; j < n; j++) {
sum -= U[i][j] * x[j];
}
x[i] = sum / U[i][i];
}
return x;
}
/**
* 交换矩阵的两行
*/
private void swapRows(double[][] matrix, int row1, int row2) {
double[] temp = matrix[row1];
matrix[row1] = matrix[row2];
matrix[row2] = temp;
}
/**
* 获取L矩阵的副本
*/
public double[][] getL() {
if (!decomposed) decompose();
double[][] copy = new double[n][n];
for (int i = 0; i < n; i++) {
System.arraycopy(L[i], 0, copy[i], 0, n);
}
return copy;
}
/**
* 获取U矩阵的副本
*/
public double[][] getU() {
if (!decomposed) decompose();
double[][] copy = new double[n][n];
for (int i = 0; i < n; i++) {
System.arraycopy(U[i], 0, copy[i], 0, n);
}
return copy;
}
/**
* 验证分解正确性: PA = LU
*/
public boolean verify(double[][] originalA) {
if (!decomposed) decompose();
double tolerance = 1e-10;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
// 计算PA[i][j] = A[pivot[i]][j]
double pa = originalA[pivot[i]][j];
// 计算LU[i][j]
double lu = 0;
for (int k = 0; k < n; k++) {
lu += L[i][k] * U[k][j];
}
if (Math.abs(pa - lu) > tolerance) {
return false;
}
}
}
return true;
}
}
矩阵运算工具：MatrixOps.java
/**
* 矩阵运算辅助类
*/
public class MatrixOps {
/**
* 矩阵与向量乘法
*/
public static double[] multiply(double[][] A, double[] x) {
int n = A.length;
double[] result = new double[n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
result[i] += A[i][j] * x[j];
}
}
return result;
}
/**
* 矩阵乘法
*/
public static double[][] multiply(double[][] A, double[][] B) {
int n = A.length;
double[][] C = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
return C;
}
/**
* 打印矩阵
*/
public static void print(double[][] A) {
for (double[] row : A) {
for (double val : row) {
System.out.printf("%10.4f ", val);
}
System.out.println();
}
}
/**
* 打印向量
*/
public static void print(double[] v) {
System.out.print("[");
for (int i = 0; i < v.length; i++) {
System.out.printf("%.4f", v[i]);
if (i < v.length - 1) System.out.print(", ");
}
System.out.println("]");
}
/**
* 生成随机矩阵
*/
public static double[][] randomMatrix(int n, java.util.Random rand) {
double[][] A = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
A[i][j] = rand.nextDouble() * 10 - 5;
}
}
// 对角占优以确保可逆
for (int i = 0; i < n; i++) {
A[i][i] = Math.abs(A[i][i]) + n * 5;
}
return A;
}
/**
* 生成Hilbert矩阵（著名的病态矩阵）
*/
public static double[][] hilbertMatrix(int n) {
double[][] H = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
H[i][j] = 1.0 / (i + j + 1);
}
}
return H;
}
}
测试入口：Main.java
public class Main {
public static void main(String[] args) {
System.out.println("=== LU分解线性代数工具测试 ===\n");
// 测试1：基本LU分解与方程组求解
System.out.println("--- 测试1：3x3方程组求解 ---");
double[][] A1 = {
{2, 1, -1},
{-3, -1, 2},
{-2, 1, 2}
};
double[] b1 = {8, -11, -3};
System.out.println("矩阵 A:");
MatrixOps.print(A1);
System.out.println("右端项 b:");
MatrixOps.print(b1);
LUDecomposition lu1 = new LUDecomposition(A1);
lu1.decompose();
System.out.println("L矩阵:");
MatrixOps.print(lu1.getL());
System.out.println("U矩阵:");
MatrixOps.print(lu1.getU());
System.out.println("分解验证: " + lu1.verify(A1));
double[] x1 = lu1.solve(b1);
System.out.println("解 x:");
MatrixOps.print(x1);
// 验证 Ax = b
double[] check = MatrixOps.multiply(A1, x1);
System.out.println("验证 Ax:");
MatrixOps.print(check);
// 测试2：行列式计算
System.out.println("\n--- 测试2：行列式计算 ---");
double[][] A2 = {
{4, 3, 2},
{1, 2, 3},
{3, 1, 4}
};
LUDecomposition lu2 = new LUDecomposition(A2);
System.out.println("矩阵:");
MatrixOps.print(A2);
System.out.printf("行列式: %.4f%n", lu2.determinant());
// 测试3：矩阵求逆
System.out.println("\n--- 测试3：矩阵求逆 ---");
double[][] A3 = {
{4, 7},
{2, 6}
};
LUDecomposition lu3 = new LUDecomposition(A3);
double[][] inv3 = lu3.inverse();
System.out.println("原矩阵:");
MatrixOps.print(A3);
System.out.println("逆矩阵:");
MatrixOps.print(inv3);
System.out.println("验证 A * A^(-1):");
MatrixOps.print(MatrixOps.multiply(A3, inv3));
// 测试4：多右端项求解
System.out.println("\n--- 测试4：多右端项求解 ---");
double[][] B = {
{8, 4},
{-11, -5},
{-3, -1}
};
double[][] X = lu1.solveMultiple(B);
System.out.println("右端项矩阵 B:");
MatrixOps.print(B);
System.out.println("解矩阵 X:");
MatrixOps.print(X);
// 测试5：大规模性能测试
System.out.println("\n--- 测试5：性能测试 (500x500) ---");
java.util.Random rand = new java.util.Random(42);
int n = 500;
double[][] bigA = MatrixOps.randomMatrix(n, rand);
double[] bigB = new double[n];
for (int i = 0; i < n; i++) bigB[i] = rand.nextDouble() * 10;
long start = System.currentTimeMillis();
LUDecomposition bigLU = new LUDecomposition(bigA);
bigLU.decompose();
long decompTime = System.currentTimeMillis() - start;
System.out.println("LU分解耗时: " + decompTime + "ms");
start = System.currentTimeMillis();
for (int i = 0; i < 100; i++) {
bigLU.solve(bigB);
}
long solveTime = System.currentTimeMillis() - start;
System.out.println("100次求解耗时: " + solveTime + "ms (平均" + (solveTime / 100.0) + "ms/次)");
// 测试6：病态矩阵
System.out.println("\n--- 测试6：Hilbert矩阵（病态）---");
double[][] H = MatrixOps.hilbertMatrix(5);
double[] Hb = {1, 0, 0, 0, 0};
System.out.println("5阶Hilbert矩阵:");
MatrixOps.print(H);
LUDecomposition luH = new LUDecomposition(H);
double[] Hx = luH.solve(Hb);
System.out.println("解:");
MatrixOps.print(Hx);
System.out.printf("行列式: %.10f%n", luH.determinant());
}
}
复杂度分析
| 操作 | 时间复杂度 | 说明 |
|---|---|---|
| LU分解 | (O(n^3)) | 一次性开销 |
| 前向替换 | (O(n^2)) | 解下三角 |
| 后向替换 | (O(n^2)) | 解上三角 |
| 单次方程组求解 | (O(n^2)) | 分解后 |
| 矩阵求逆 | (O(n^3)) | n次求解 |
| 行列式计算 | (O(n)) | 分解后U对角线乘积 |
空间复杂度：(O(n^2))，存储L和U矩阵。
数值稳定性
部分主元选取（Partial Pivoting）是LU分解数值稳定的关键：
– 每次消元前在下方选取绝对值最大的元素作为主元
– 将乘数 (m_{ik}) 的绝对值控制在 (\leq 1)
– 有效抑制舍入误差的传播和放大
对于绝大多数实际矩阵，部分主元LU分解的数值稳定性足以满足工程需求。仅在极少数特殊病态矩阵（如Hilbert矩阵）上需要更复杂的完全主元选取。
适用场景
- 工程计算：结构力学有限元分析、电路网络分析
- 优化算法：线性规划、最小二乘问题的正规方程求解
- 计算机图形学：变换矩阵求逆、投影计算
- 机器学习：岭回归、高斯过程协方差矩阵求解
- 经济学：投入产出分析、一般均衡模型
与QR分解、Cholesky分解的对比
| 分解方法 | 适用矩阵 | 稳定性 | 计算量 |
|---|---|---|---|
| LU | 一般方阵 | 中（需主元） | (O(n^3/3)) |
| QR | 一般矩阵 | 高 | (O(2n^3/3)) |
| Cholesky | 对称正定 | 高 | (O(n^3/6)) |
扩展方向
- 块LU分解：利用缓存局部性加速大规模矩阵
- 稀疏LU：针对稀疏矩阵的符号分解与数值分解分离
- 不完全LU（ILU）：预处理共轭梯度法的重要工具
- PLU分解：将置换矩阵P显式维护，优化批量求解
总结
LU分解是数值线性代数的基石算法。通过Java实现，我们理解了将矩阵分解为下三角和上三角因子的核心思想，以及部分主元选取对数值稳定性的关键作用。在需要多次求解相同系数矩阵的不同右端项时，LU分解通过一次 (O(n^3)) 分解加多次 (O(n^2)) 求解，大幅提升了计算效率。从工程有限元到机器学习，从图形学到经济学，LU分解是每一位算法工程师必须掌握的基础工具。
