Gauss-Jordan 消元法复习笔记
因为洛谷全站推荐审核问题，洛谷专栏可能无法及时更新。
Gauss-Jordan 消元法用于解决多元一次方程组求根问题。
名词解释
-
高斯消元（Gaussian Elimination）：用于求解多元一次方程组求根问题
-
Gauss-Jordan 消元（高斯-约旦消元，高斯-乔丹消元）：高斯消元的一种，但是代码更好写，且可以一次性求出最简行阶梯形矩阵/规范行阶梯形矩阵
构建增广矩阵
在高斯消元中，一般会把 \(n\) 个 \(n\) 元一次方程组写成一个 \(n\) 行 \((n+1)\) 列的“增广矩阵”的形式，具体来说：
- 把方程移项以格式化为：左边全部都是多个“系数 \(\times x_i\)”的和的形式，若某一个方程中 \(x_i\) 不存在则给他补上一个 \(0\) 的系数；常数项则全部移动到右侧。
- 这样，每一个方程等号左侧恰有 \(n\) 个一次项（系数为 \(0\) 也算作一次项），右侧恰有一个常数项。
- 矩阵的第 \(i\) 行 \(j\) 列表示第 \(i\) 个方程中 \(x_j\) 的系数；特别的，第 \(i\) 行 \((n+1)\) 列表示方程右侧的常数项。
（其实通常情况下题目中会直接把增广矩阵给你，而不需要手动处理。）
举例：
第一步：方程格式化（移项，补 \(0\) 系数）
第二步：构建增广矩阵（直接把上面的系数做成一个矩阵）
竖线左边表示系数，右边表示常数。此后高斯消元就是在这个矩阵上进行操作。
消元
我们可以对这个矩阵进行如下操作：
- 对不同行调换位置
- 对某一行全部乘以或除以一个常数
- 让某一行所有数加上或减去另一行所有数
任意进行这三个操作后的矩阵与原矩阵等价，并且只用这三个操作就可以完成整个高斯消元过程。
Gauss-Jordan 消元法的特点是可以把一个合法的方程组一次性消成“单位矩阵”（除了竖线右边的增广部分，左侧的主对角线全部都是 \(1\)）。
具体来说，假设矩阵第 \(x\) 行第 \(y\) 列为 \(M(x, y)\)，步骤如下：
-
按照列依次处理，对于第 \(i\) 列，找到某个 \(p \ge i\)，使 \(M(p, i)\) 不为 \(0\)。
-
把 \(p\) 行与第 \(i\) 交换（可以大大简化后续代码书写）。于是我们可以只用这一行来将其它所有行的第 \(i\) 列都变为 \(0\)。
-
先让第 \(i\) 行所有数都除以 \(M(i, i)\)，这样可以让所有 \(M(i, i)\) 变为 \(1\) 且矩阵与原矩阵等价。
-
枚举所有其它行，即对于每一个 \(j \neq i\)，把这一整行都减去 \(M(j, i)\) 乘以第 \(i\) 行，于是 \(M(j, i)\) 就被消去了。
（其实应该是用 \(M(j, i) \div M(i, i)\) 来乘，不过 \(M(i, i) = 1\) 被约掉了） -
这样，第 \(i\) 列除了第 \(i\) 行之外，全部被清零。对于每一个 \(i \in [1, n]\) 都做一遍，最后得到的是一个“最简行阶梯形矩阵”或“规范行阶梯形矩阵”，即矩阵竖线左边就是一个单位矩阵，其增广部分直接就是方程组的根了。
继续使用上面的示例。
- \(\textcolor{Blue}{蓝色}\)表示已经处理完毕保证不会改变的数
- \(\textcolor{Green}{绿色}\)表示处理完毕且保证会是 \(0\) 的数
- \(\textcolor{Orange}{橙色}\)表示正在处理的列
- \(\textcolor{Gold}{黄色}\)表示用来消去其它行的行
- \(\textcolor{Red}{红色}\)表示黄色和橙色的相交处——同时也是这一次消元的主角
第一列
（第 \(1\)、\(2\) 步）首先是第 \(1\) 列，我们发现第 \(1\) 行就不是 \(0\)，可以拿来消元。于是把它换到第 \(1\) 行上来（好吧，其实相当于啥也没干）。
（第 \(3\) 步）然后让第 \(1\) 行所有数全部除以 \(\textcolor{red}{3}\)，目的是把 \(\textcolor{red}{3}\) 变成 \(\textcolor{blue}{1}\)（它以后不会再改变了）。
（第 \(4\) 步）第 \(2\) 行减去 \(\textcolor{Orange}{1}\) 倍的第 \(1\) 行，第 \(3\) 行减去 \(\textcolor{Orange}{4}\) 倍的第一行，这样这两者橙色数字都一定会被消成 \(\textcolor{Green}{0}\)。
第二列
（第 \(1\)、\(2\) 步）现在要来处理第 \(2\) 列。我们找到 \(M(3, 2)\) 是不为 \(0\) 的（也是唯一一个合法的），于是把第 \(3\) 行换上来：
（第 \(3\) 步）第 \(2\) 行全部除以 \(\textcolor{Red}{2}\)。
（第 \(4\) 步）第 \(1\) 行减去 \(\textcolor{Orange}{0}\) 倍的第 \(2\) 行，第 \(3\) 行减去 \(\textcolor{Orange}{0}\) 倍的第 \(2\) 行（其实都相当于啥也没干，不过形式还是要走一下的），确保这两个橙色数字被消成 \(\textcolor{Green}{0}\)。
第三列
（第 \(1\)、\(2\) 步）第 \(3\) 列，我们找到 \(M(3, 3)\) 是不为 \(0\) 的（也是唯一一个合法的，因为其它的要么为 \(0\)，要么在已经处理一次过的行里面），所以把它换到第 \(3\) 行（其实相当于啥也没干）。
（第 \(3\) 步）第 \(3\) 行全部除以 \(\textcolor{Red}{\frac{11}{3}}\)。
（第 \(4\) 步）第 \(1\) 行减去 \(\textcolor{Orange}{\frac{4}{3}}\) 倍的第 \(3\) 行，第 \(2\) 行减去 \(\textcolor{Orange}{-\frac{7}{6}}\) 倍的第 \(3\) 行，使得这两个橙色数字全部变成 \(\textcolor{Green}{0}\)。
最后
可以看到，这个矩阵竖线左边的部分已经是一个单位矩阵了，我们把它还原成方程：
简化：
这就是方程的根。
还有一些细节：
为什么必须保证 \(p \ge i\)？
因为在只有 \(i\) 及以后的行才能保证这一行的第 \(1\) 到 \((i - 1)\) 列全都为 \(0\)，和前面的行做加减时不会导致已经形成的部分单位矩阵再次被破坏。
要是没有找到合法的 \(p\) 呢？
这说明原方程组没有唯一合法解，具体在接下来会讲。
无解处理
上述情况都是有解时高斯消元的正常处理过程，那么怎么判断无解呢？
这里所说的无解其实是“没有唯一合法解”，即两种情况：“没有合法解”和“有无穷多合法解”。
这里需要一些线性代数知识（“线性相关”知识）。不过简单来说就是：如果在某一列发现没有合法的 \(p\)，即第 \(i\) 行及后面全都是 \(0\)，那就说明这里面有一个式子是可以通过其它式子凑出来的——它不是一个真正有用的式子。此时相当于 \(n\) 个未知数但是只有 \((n-1)\) 乃至更少的方程式，方程自然就无唯一合法解了。
此时可以直接判定“没有唯一合法解”了。但是有些题目会要求判断到底是“没有合法解”还是“有无穷多合法解”，这就需要对逻辑进行一定改动了。
具体来说，我们并不追求“完美消元”，而是追求“尽量消元”，当这一列全部找不到合法 \(p\) 的时候，我们就直接跳过这一列，在下一列继续匹配。
这样，最后若是没有“完美匹配”，有些列无法找到第 \(p\) 行匹配，一定会剩下一些行堆积在矩阵下方，且这些行的系数可以保证全部为 \(0\)（对于这一行的每一列，如果这一列曾匹配成功会被消成 \(0\)，若没有匹配成功说明找不到位置合法且非 \(0\) 的 \(p\)，这里一定还是 \(0\)）。
对于这些式子，把他们单独提取出来：
如果所有上面出现的所有 \(v\) 都是 \(0\)，那么前面的 \(x\) 无论取什么都可以——有无穷多解；如果上面的 \(v\) 中出现了非 \(0\) 数，那么 \(x\) 怎么取都没法凑出 \(0\) 以外的数字来——无合法解。
另外，如果题目是给了你超过 \(n\) 个方程，也可以用这种方式来判断。
代码
在代码中要注意：因为大量实数运算必定会产生精度问题，所以对 \(0\) 的判断不能用 \(v = 0\) 而应该用 \(\lvert v \rvert < \operatorname{EPS}\)，其中 \(\operatorname{EPS}\) 表示一个极小数，如 \(10^{-6}\)。
下面给出“洛谷 P2455 [SDOI2006] 线性方程组”一题的参考代码。
#include <iomanip>
#include <iostream>
using namespace std;
constexpr int MAXN = 55;
int n; double a[MAXN][MAXN];
const double EPS = 1e-6;
int Gauss_Jordan()
{
int row = 1; // 当前还未处理的最早的行数
for (int col = 1; col <= n; col++) // 枚举列数
{
/* 第 1 步 */
int p = row; // 在没处理过的行里找 p
while (p <= n && abs(a[p][col]) < EPS) p++;
if (p > n) continue; // 找不到合法 p，跳过这一行
/* 第 2 步 */
for (int i = 1; i <= n + 1; i++)
swap(a[p][i], a[row][i]); // 交换两行
/* 第 3 步 */
double tmp = a[row][col]; // 预存，因为后面会修改其值
for (int i = 1; i <= n + 1; i++)
a[row][i] /= tmp; // 把 a[row][col] 调成 1
/* 第 4 步 */
for (int i = 1; i <= n; i++)
{
if (i == row) continue; // 不能把自己消了
tmp = a[i][col]; // （同上）预存，因为后面会修改其值
for (int j = 1; j <= n + 1; j++)
a[i][j] -= a[row][j] * tmp; // 目的是消掉 a[i][row]
}
row++;
}
if (row == n + 1) return 1; // 有唯一合法解
for (int i = row; i <= n; i++)
if (abs(a[i][n + 1]) >= EPS) return -1; // 无合法解
return 0; // 有无穷多合法解
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
cin >> n;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n + 1; j++)
cin >> a[i][j];
int ans = Gauss_Jordan();
if (ans == 1)
{
for (int i = 1; i <= n; i++)
cout << 'x' << i << '=' << setprecision(3) << a[i][n + 1] << '\n';
}
else cout << ans << '\n';
return 0;
}
“洛谷 P3389 【模板】高斯消元法”的代码不再给注释，且只给出关键部分（这个代码是我去年打的，所以码风略有不同）
const double EPS=1e-6;
bool Gauss_Jordan()
{
for(int i=1;i<=n;i++)
{
int p=i;
while(abs(a[p][i])<EPS && p<=n) p++;
if(p>n) return false;
for(int j=1;j<=n+1;j++)
swap(a[i][j],a[p][j]);
double si=a[i][i];
for(int j=i;j<=n+1;j++)
a[i][j]/=si;
for(int j=1;j<=n;j++)
{
if(j==i) continue;
double sj=a[j][i];
for(int k=i;k<=n+1;k++)
a[j][k]-=sj*a[i][k];
}
}
return true;
}
本文采用 「CC-BY-NC 4.0」 创作共享协议，转载请注明作者及出处，禁止商业使用。
作者：Jerrycyx，原文链接：https://www.cnblogs.com/jerrycyx/p/19073499
浙公网安备 33010602011771号