每日算法 — 使用java实现高斯消元法：线性方程组求解与矩阵行变换

算法背景

高斯消元法（Gaussian Elimination）是线性代数中求解线性方程组、计算矩阵秩、求逆矩阵的核心算法。由Carl Friedrich Gauss系统化提出，其基本思想是通过初等行变换将增广矩阵化为行阶梯形（前向消元），再通过回代求解未知数。该方法在计算机图形学（坐标变换）、经济学（投入产出分析）、工程计算（结构力学）等领域不可或缺。

核心思想

高斯消元分为两个阶段：
1. 前向消元：选取主元，通过行变换将主元下方元素消为0，最终得到上三角矩阵
2. 回代求解：从最后一行开始，依次代入已求得的变量值，解出当前未知数

为提升数值稳定性，实践中常采用列主元消去法（每列选取绝对值最大的元素作为主元）。

Java实现

import java.util.Arrays;

/**
 * 高斯消元法实现
 * 支持线性方程组求解、矩阵求逆与秩计算
 */
public class GaussianElimination {

    private static final double EPS = 1e-9;

    /**
     * 使用高斯消元法求解线性方程组 Ax = b
     * @param a 系数矩阵（n x n）
     * @param b 常数向量（n）
     * @return 解向量x，若无解或多解返回null
     */
    public static double[] solve(double[][] a, double[] b) {
        int n = a.length;
        // 构造增广矩阵
        double[][] aug = new double[n][n + 1];
        for (int i = 0; i < n; i++) {
            System.arraycopy(a[i], 0, aug[i], 0, n);
            aug[i][n] = b[i];
        }

        // 前向消元
        for (int col = 0, row = 0; col < n && row < n; col++, row++) {
            // 列主元选取：在当前列下方找绝对值最大的元素
            int pivotRow = row;
            for (int i = row + 1; i < n; i++) {
                if (Math.abs(aug[i][col]) > Math.abs(aug[pivotRow][col])) {
                    pivotRow = i;
                }
            }

            // 如果主元近似为0，说明矩阵奇异
            if (Math.abs(aug[pivotRow][col]) < EPS) {
                return null; // 无解或无穷多解
            }

            // 交换当前行与主元行
            double[] temp = aug[row];
            aug[row] = aug[pivotRow];
            aug[pivotRow] = temp;

            // 将主元归一化为1（可选，不影响结果但便于观察）
            double pivot = aug[row][col];
            for (int j = col; j <= n; j++) {
                aug[row][j] /= pivot;
            }

            // 消去下方所有行的当前列元素
            for (int i = row + 1; i < n; i++) {
                double factor = aug[i][col];
                if (Math.abs(factor) < EPS) continue;
                for (int j = col; j <= n; j++) {
                    aug[i][j] -= factor * aug[row][j];
                }
            }
        }

        // 回代求解
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; i--) {
            x[i] = aug[i][n];
            for (int j = i + 1; j < n; j++) {
                x[i] -= aug[i][j] * x[j];
            }
            // 如果系数矩阵已归一化，此处无需再除
            if (Math.abs(aug[i][i]) > EPS) {
                x[i] /= aug[i][i];
            }
        }

        return x;
    }

    /**
     * 计算矩阵的秩
     */
    public static int rank(double[][] matrix) {
        int rows = matrix.length;
        int cols = matrix[0].length;
        double[][] m = new double[rows][cols];
        for (int i = 0; i < rows; i++) {
            m[i] = matrix[i].clone();
        }

        int r = 0;
        for (int c = 0; c < cols && r < rows; c++) {
            // 找主元
            int pivotRow = r;
            for (int i = r + 1; i < rows; i++) {
                if (Math.abs(m[i][c]) > Math.abs(m[pivotRow][c])) {
                    pivotRow = i;
                }
            }

            if (Math.abs(m[pivotRow][c]) < EPS) continue;

            // 交换
            double[] temp = m[r];
            m[r] = m[pivotRow];
            m[pivotRow] = temp;

            // 消元
            for (int i = r + 1; i < rows; i++) {
                double factor = m[i][c] / m[r][c];
                for (int j = c; j < cols; j++) {
                    m[i][j] -= factor * m[r][j];
                }
            }
            r++;
        }

        return r;
    }

    /**
     * 计算矩阵的行列式（仅方阵）
     */
    public static double determinant(double[][] matrix) {
        int n = matrix.length;
        double[][] m = new double[n][n];
        for (int i = 0; i < n; i++) {
            m[i] = matrix[i].clone();
        }

        double det = 1.0;
        int swapCount = 0;

        for (int i = 0; i < n; i++) {
            // 找主元
            int pivotRow = i;
            for (int j = i + 1; j < n; j++) {
                if (Math.abs(m[j][i]) > Math.abs(m[pivotRow][i])) {
                    pivotRow = j;
                }
            }

            if (Math.abs(m[pivotRow][i]) < EPS) {
                return 0; // 奇异矩阵，行列式为0
            }

            if (pivotRow != i) {
                double[] temp = m[i];
                m[i] = m[pivotRow];
                m[pivotRow] = temp;
                swapCount++;
            }

            det *= m[i][i];

            // 消去下方元素
            for (int j = i + 1; j < n; j++) {
                double factor = m[j][i] / m[i][i];
                for (int k = i; k < n; k++) {
                    m[j][k] -= factor * m[i][k];
                }
            }
        }

        // 行交换改变行列式符号
        if (swapCount % 2 == 1) {
            det = -det;
        }

        return det;
    }

    /**
     * 矩阵求逆（使用增广矩阵 [A|I] 变换为 [I|A^-1]）
     */
    public static double[][] inverse(double[][] matrix) {
        int n = matrix.length;

        // 构造增广矩阵 [A | I]
        double[][] aug = new double[n][2 * n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(matrix[i], 0, aug[i], 0, n);
            aug[i][n + i] = 1.0;
        }

        // 高斯-约旦消元
        for (int i = 0; i < n; i++) {
            // 列主元
            int pivotRow = i;
            for (int j = i + 1; j < n; j++) {
                if (Math.abs(aug[j][i]) > Math.abs(aug[pivotRow][i])) {
                    pivotRow = j;
                }
            }

            if (Math.abs(aug[pivotRow][i]) < EPS) {
                return null; // 矩阵不可逆
            }

            // 交换
            double[] temp = aug[i];
            aug[i] = aug[pivotRow];
            aug[pivotRow] = temp;

            // 归一化当前行
            double pivot = aug[i][i];
            for (int j = 0; j < 2 * n; j++) {
                aug[i][j] /= pivot;
            }

            // 消去其他行的当前列
            for (int j = 0; j < n; j++) {
                if (j == i) continue;
                double factor = aug[j][i];
                if (Math.abs(factor) < EPS) continue;
                for (int k = 0; k < 2 * n; k++) {
                    aug[j][k] -= factor * aug[i][k];
                }
            }
        }

        // 提取逆矩阵
        double[][] inv = new double[n][n];
        for (int i = 0; i < n; i++) {
            System.arraycopy(aug[i], n, inv[i], 0, n);
        }

        return inv;
    }

    /**
     * 打印矩阵
     */
    public static void printMatrix(double[][] matrix) {
        for (double[] row : matrix) {
            for (double val : row) {
                System.out.printf("%10.4f ", val);
            }
            System.out.println();
        }
    }

    /**
     * 打印向量
     */
    public static void printVector(double[] vec) {
        System.out.println(Arrays.toString(vec));
    }

    public static void main(String[] args) {
        System.out.println("=== 线性方程组求解 ===");
        // 方程组：
        // 2x + y - z = 8
        // -3x - y + 2z = -11
        // -2x + y + 2z = -3
        double[][] a = {
            {2, 1, -1},
            {-3, -1, 2},
            {-2, 1, 2}
        };
        double[] b = {8, -11, -3};

        System.out.println("系数矩阵 A:");
        printMatrix(a);
        System.out.println("常数向量 b: " + Arrays.toString(b));

        double[] x = solve(a, b);
        if (x != null) {
            System.out.println("解向量 x: " + Arrays.toString(x));
            // 验证
            System.out.println("验证 Ax:");
            for (int i = 0; i < a.length; i++) {
                double sum = 0;
                for (int j = 0; j < a[0].length; j++) {
                    sum += a[i][j] * x[j];
                }
                System.out.printf("  方程 %d: %.4f (应为 %.4f)%n", i + 1, sum, b[i]);
            }
        } else {
            System.out.println("方程组无解或有无穷多解");
        }

        System.out.println("\n=== 矩阵求逆 ===");
        double[][] inv = inverse(a);
        if (inv != null) {
            System.out.println("逆矩阵 A^-1:");
            printMatrix(inv);
        } else {
            System.out.println("矩阵不可逆");
        }

        System.out.println("\n=== 行列式计算 ===");
        double det = determinant(a);
        System.out.printf("det(A) = %.4f%n", det);

        System.out.println("\n=== 矩阵秩计算 ===");
        double[][] rankMatrix = {
            {1, 2, 3},
            {2, 4, 6},
            {1, 1, 1}
        };
        System.out.println("矩阵:");
        printMatrix(rankMatrix);
        System.out.println("秩: " + rank(rankMatrix));

        System.out.println("\n=== 奇异矩阵检测 ===");
        double[][] singular = {
            {1, 2, 3},
            {4, 5, 6},
            {7, 8, 9}
        };
        System.out.println("矩阵:");
        printMatrix(singular);
        System.out.printf("行列式: %.4f%n", determinant(singular));
        System.out.println("逆矩阵: " + (inverse(singular) == null ? "不存在" : "存在"));
    }
}

复杂度分析

算法优势与适用场景

优势：
– 通用性强，可处理任意规模的稠密线性方程组
– 列主元消去大幅提升数值稳定性
– 通过增广矩阵扩展可统一处理求逆、求秩等问题

适用场景：
– 工程计算中的结构力学与电路分析
– 计算机图形学的三维坐标变换
– 机器学习中的最小二乘线性回归
– 经济学投入产出模型求解