高斯消元学习笔记
引入
高斯-约旦消元法（Gauss–Jordan elimination）是求解线性方程组的经典算法，它在当代数学中有着重要的地位和价值，是线性代数课程教学的重要组成部分。
高斯消元法除了用于线性方程组求解外，还可以用于行列式计算、求矩阵的逆，以及其他计算机和工程方面。
过程
一个经典的问题，给定一个有 \(n\) 个未知数的线性方程组，求方程组的解。
例如：
因为未知数的值在运算中也不会发生变化，可以尝试用矩阵将其转换。
定义
增广矩阵
一个矩阵，如果在矩阵右方有一个列向量或者一个列向量矩阵，那么这个矩阵就是增广矩阵。
将上方的例子改写成增广矩阵就是下图：
可以发现最终的答案需要一个对角矩阵。
那么怎么将增广矩阵变换成对角矩阵呢。
将增广矩阵变成对角矩阵的过程
假定当前求解第 \(i\) 个未知数 \(x_i\)。
过程如下：
- 找到第一行 \(j\) 满足矩阵 \(mat\) 的 \(mat_{j,i} \neq 0 \land (mat_{j,j} \neq 1 \lor j \ge i)\) 。
- 交换第 \(i\) 行和第 \(j\) 行。
- 将其他行的第 \(i\) 列变成 \(0\)。
- 当 \(i \le n\) 时，重复以上步骤。
矩阵的秩
一个经过高斯消元的对角矩阵 \(A\) 如果满足第 \(i\) 行不是全零且第 \(i + 1\) 为全零，则矩阵的秩为 \(i\)，记作 \(r(A)\)。
无解情况
因为最后的结果会是一个对角方阵且此方阵的秩为 \(n\) 且满足此矩阵 \(A\) 中 \(\exists i, A_{i,i} = 0\) 并且右方常数不为 \(0\)。
无数解情况
满足一个对角矩阵 \(A\)，\(r(A) = n \land \exists i,A_{i,i} = 0\) 且其对应的常数为 \(0\)。
可以发现如果一个矩阵有唯一解满足 \(r(A) = n \land \forall i A_{i,i} \ne 0\)。
实现
一下代码为 SDOI 2006 线性方程组 的代码。
#include <cstdio>
#include <cmath>
#include <algorithm>
using u32 = unsigned int ;
using i64 = long long ;
using u64 = unsigned long long ;
using f32 = float ;
using f64 = double ;
const f64 eps = 1e-6 ;
const int N = 105 ;
int n;
f64 mat[N][N];
void Gauss(){
for (int i = 1; i <= n; i++) {
int p = 0;
for (int j = 1; j <= n; j++) {
if (fabs(mat[j][j]) > eps && j < i) {
continue;
}
if (fabs(mat[j][i]) > fabs(mat[p][i])) {
p = j;
break;
}
}
if (!p) {
continue;
}
std::swap(mat[i], mat[p]);
for (int j = 1; j <= n; j++) {
if (i == j) {
continue;
}
double val = mat[j][i] / mat[i][i];
for (int k = 1; k <= n + 1; k++) {
mat[j][k] -= mat[i][k] * val;
}
}
}
}
int main(){
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n + 1; j++) {
scanf("%lf", &mat[i][j]);
}
}
Gauss();
for (int i = 1; i <= n; i++) {
if (fabs(mat[i][n + 1]) > eps && fabs(mat[i][i]) < eps) {
puts("-1");
return 0;
}
}
for (int i = 1; i <= n; i++) {
if (fabs(mat[i][n + 1]) < eps && fabs(mat[i][i]) < eps) {
puts("0");
return 0;
}
}
for (int i = 1; i <= n; i++) {
double ans = mat[i][n + 1] / mat[i][i];
if (fabs(ans) < eps) {
puts("0.00");
} else {
printf("%.2lf\n", ans);
}
}
return 0;
}
时间复杂度
显然，\(\mathcal{O}(n ^ 3)\)。
其他应用
行列式求值
一个行列式如果它是对角矩阵，那么它的值为对角的乘积。
而高斯消元正好能将一个矩阵变换成对角矩阵，因此可以在 \(\mathcal{O}(n ^ 3)\) 求解。
实现
以下为 【模板】行列式求值 代码，通过辗转相减法实现了求值取模 \(p\) 的答案。
#include <cstdio>
#include <algorithm>
using i64 = long long ;
const int N = 605 ;
int n, mod;
i64 a[N][N];
i64 ans = 1;
void Gauss(){
for (int i = 1; i <= n; i++) {
for (int j = i + 1; j <= n; j++) {
while (a[i][i]) {
int t = a[j][i] / a[i][i] ;
for (int k = 1; k <= n; k++) {
a[j][k] = (a[j][k] - t * a[i][k] % mod + mod) % mod;
}
ans *= -1;//行列式的性质，交换两行，需要改变结果的符号
std::swap(a[i], a[j]);
}
ans *= - 1;
std::swap(a[i], a[j]);
}
}
}
int main(){
scanf("%d %d", &n, &mod);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
scanf("%lld", &a[i][j]);
}
}
Gauss();
for (int i = 1; i <= n; i++) {
ans = (ans * a[i][i] + mod) % mod;
}
printf("%lld",(ans + mod) % mod);
return 0;
}
矩阵求逆
矩阵求逆就是通过高斯消元将此矩阵消成单位矩阵，而右方可以记录矩阵的初等变换的积，最后右方就是矩阵的逆。
实现
以下为 【模板】矩阵求逆 的代码。
#include <cstdio>
#include <cstring>
#include <algorithm>
using u32 = unsigned int;
using i64 = long long ;
using u64 = unsigned long long ;
const int mod = 1e9 + 7 ;
const int N = 405 ;
int n;
int mat[N][N << 1];
i64 qpow(i64 a, i64 b = mod - 2) {
i64 ans = 1;
for (; b; b >>= 1, a = a * a % mod) {
if (b & 1) {
ans = ans * a % mod;
}
}
return ans;
}
int Gauss(){
for (int i = 1; i <= n; i++) {
int p = 0;
for (int j = 1; j <= n; j++) {
if (j < i && mat[j][j]) {
continue;
}
if (mat[j][i]) {
p = j;
break;
}
}
if (!p) {
return -1;//如果方程组无解，自然就没有
}
std::swap(mat[i], mat[p]);
i64 inv = qpow(mat[i][i]);
for (int j = 1; j <= n; j++) {
if (i == j) {
continue;
}
i64 val = mat[j][i] * inv % mod;
for (int k = 1; k <= (n << 1); k++) {
mat[j][k] = (mat[j][k] - val * mat[i][k] % mod + mod) % mod;
}
}
for (int j = 1; j <= (n << 1); j++) {
mat[i][j] = (mat[i][j] * inv) % mod;
}
}
return 1;
}
int main(){
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
scanf("%d", &mat[i][j]);
mat[i][n + i] = 1;
}
}
if (Gauss() == -1) {
printf("No Solution");
return 0;
}
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
printf("%d ", mat[i][n + j]);
}
puts("");
}
return 0;
}
求解模意义下的方程组
很简单，将除法改成逆元就行了。
特殊的，二进制下可以用异或操作代替逆元。
