1. 引言
高斯消元法（Gaussian Elimination）是线性代数中最基础、最重要的算法之一，用于求解线性方程组、计算矩阵的秩、求逆矩阵以及计算行列式等。该方法由德国数学家卡尔·弗里德里希·高斯在 19 世纪初系统化，但类似的思想早在古代中国《九章算术》中就已出现。
高斯消元法的核心思想是通过一系列初等行变换，将系数矩阵化为行阶梯形（Row Echelon Form）或简化行阶梯形（Reduced Row Echelon Form），从而简化方程组的求解过程。
2. 算法原理
2.1 线性方程组表示
一个包含 n 个未知数、m 个方程的线性方程组可以表示为：
a₁₁x₁ + a₁₂x₂ + … + a₁ₙxₙ = b₁
a₂₁x₁ + a₂₂x₂ + … + a₂ₙxₙ = b₂
…
aₘ₁x₁ + aₘ₂x₂ + … + aₘₙxₙ = bₘ
用矩阵形式表示为：Ax = b，其中：
- A 是 m×n 的系数矩阵；
- x 是 n×1 的未知数向量；
- b 是 m×1 的常数向量。
2.2 增广矩阵
为了方便计算，我们将系数矩阵 A 和常数向量 b 合并为增广矩阵（augmented matrix）：
[A | b] =
[ a₁₁ a₁₂ … a₁ₙ | b₁ ]
[ a₂₁ a₂₂ … a₂ₙ | b₂ ]
[ … … … … | … ]
[ aₘ₁ aₘ₂ … aₘₙ | bₘ ]
2.3 初等行变换
高斯消元法允许的三种初等行变换为：
这些变换不会改变线性方程组的解集。
3. 算法步骤
3.1 前向消元（Forward Elimination）
目标：将增广矩阵化为行阶梯形（Row Echelon Form, REF）
行阶梯形的特征：
算法步骤：
3.2 回代求解（Back Substitution）
目标：将行阶梯形化为简化行阶梯形（RREF）并求解。
简化行阶梯形的特征：
算法步骤：
4. 代码实现
4.1 Python 实现
import numpy as np
def gaussian_elimination(A, b):
"""
高斯消元法求解线性方程组 Ax = b
参数：
A: n×n 系数矩阵
b: n×1 常数向量
返回：
x: 解向量
"""
n = len(A)
# 构造增广矩阵
Ab = np.hstack([A, b.reshape(–1, 1)])
# 前向消元
for i in range(n):
# 选主元（部分选主元法）
max_row = i + np.argmax(np.abs(Ab[i:, i]))
if max_row != i:
Ab[[i, max_row]] = Ab[[max_row, i]]
# 如果主元为 0，方程组无解或有无穷多解
if np.abs(Ab[i, i]) < 1e-10:
raise ValueError("矩阵奇异或接近奇异，无法求解")
# 消去下方行的当前列元素
for j in range(i + 1, n):
factor = Ab[j, i] / Ab[i, i]
Ab[j, i:] -= factor * Ab[i, i:]
# 回代求解
x = np.zeros(n)
for i in range(n – 1, –1, –1):
x[i] = (Ab[i, –1] – np.dot(Ab[i, i+1:n], x[i+1:])) / Ab[i, i]
return x
# 示例
A = np.array([[2, 1, –1],
[–3, –1, 2],
[–2, 1, 2]], dtype=float)
b = np.array([8, –11, –3], dtype=float)
x = gaussian_elimination(A, b)
print("解向量 x =", x)
print("验证 Ax – b =", np.dot(A, x) – b)
4.2 C++ 实现
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
vector<double> gaussianElimination(vector<vector<double>> A, vector<double> b) {
int n = A.size();
// 构造增广矩阵
vector<vector<double>> Ab(n, vector<double>(n + 1));
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
Ab[i][j] = A[i][j];
}
Ab[i][n] = b[i];
}
// 前向消元
for (int i = 0; i < n; i++) {
// 选主元
int maxRow = i;
double maxVal = fabs(Ab[i][i]);
for (int k = i + 1; k < n; k++) {
if (fabs(Ab[k][i]) > maxVal) {
maxVal = fabs(Ab[k][i]);
maxRow = k;
}
}
// 交换行
if (maxRow != i) {
swap(Ab[i], Ab[maxRow]);
}
// 检查主元是否为0
if (fabs(Ab[i][i]) < 1e-10) {
throw runtime_error("矩阵奇异，无法求解");
}
// 消去下方行的当前列元素
for (int k = i + 1; k < n; k++) {
double factor = Ab[k][i] / Ab[i][i];
for (int j = i; j <= n; j++) {
Ab[k][j] -= factor * Ab[i][j];
}
}
}
// 回代求解
vector<double> x(n);
for (int i = n – 1; i >= 0; i—) {
double sum = 0.0;
for (int j = i + 1; j < n; j++) {
sum += Ab[i][j] * x[j];
}
x[i] = (Ab[i][n] – sum) / Ab[i][i];
}
return x;
}
int main() {
vector<vector<double>> A = {{2, 1, –1},
{–3, –1, 2},
{–2, 1, 2}};
vector<double> b = {8, –11, –3};
vector<double> x = gaussianElimination(A, b);
cout << "解向量 x = [";
for (int i = 0; i < x.size(); i++) {
cout << x[i];
if (i < x.size() – 1) cout << ", ";
}
cout << "]" << endl;
return 0;
}
5. 算法复杂度分析
5.1 时间复杂度
对于 n×n 的矩阵：
- 前向消元：约需要 2n³/3 次浮点运算。
- 回代求解：约需要 n² 次浮点运算。
- 总复杂度：O(n³)。
对于 m×n 的矩阵（m ≥ n）：
- 总复杂度为 O(mn²)
5.2 空间复杂度
- 原地算法：O(1) 额外空间（直接在原矩阵上操作）
- 非原地算法：O(mn) 额外空间
