【算法笔记】Gauss 消元 · 行列式
- 本文总计约 25000 字，阅读大约需要 75 分钟。
- 警告！警告！警告！本文有大量的 \(\LaTeX\) 公式渲染，可能会导致加载异常缓慢！
前言
啊啊，终于搞完了啊~Gauss 消元法相关的知识。
实际上，Gauss 消元的内容并不复杂，但是相比较于 Gauss 消元本身，其背后的线性代数的知识更加重要。所以，我决定以介绍 Gauss 消元法为契机，更多地讲一讲线性代数相关的知识。
这并不只是为了学习 OI 知识，线性代数，也是大学里必修的一门课程。
你，准备好了吗？
问题引入
『鸡兔同笼』问题是我国古代的一道经典数学问题，出自于《孙子算经》：
今有鸡兔同笼，上有三十五头，下有九十四足，问鸡兔各几何？
用现代的数学语言来表示，设鸡的数目为 \(x\)，兔的数目为 \(y\)，那么可以列如下的方程组：
这是一个二元一次方程组，当然可以求解出来 \(x=23\)，\(y=12\)。
那么，如果给定一个 \(n\) 元一次方程组，那么能否也求出它的解呢？
前置知识
矩阵
矩阵是用数字写成的，有一定行数和列数的数表。常用大写字母 \(A,B,\cdots,C\) 表示。例如，下面的这个，就是一个 \(2\) 行 \(3\) 列的矩阵。
『真是废话的定义呢……那这个矩阵，有什么用处嘛？』
呃呃，矩阵的定义是挺简单粗暴的，但是它可以启发我们求 \(n\) 元一次方程组的解。
考虑下面的这个方程组：
我们可以用矩阵的方式来表示这个方程组，如下：
像这样的矩阵，我们就将其称为一个方程组的增广矩阵。
初等行变换
对于一个矩阵，定义以下的三个操作为它的初等行变换：
- 交换某两行，例如：
- 将某一行内所有元素均乘以一个非 \(0\) 的常数 \(k\)，例如：
- 将某一行与另一行相加，例如：
『那么，如果对一个方程组的增广矩阵做初等行变换的话……
就用上面「鸡兔同笼」的方程组为例吧，它的增广矩阵就是：
如果我们把这个矩阵的第一行加到第二行上的话，就变成了这样：
诶诶，这两个矩阵对应的方程是同解的欸？
那么是不是任何一个增广矩阵，在做初等行变换以后，得到的方程与原方程都是同解的呢？』
答案是肯定的！这个命题很显然，做初等行变换并不会改变方程组的解。所以，这个优秀的性质也为我们介绍 Gauss 消元法奠定了基础。
不过在此之前，我们还得再解释一些定义。
矩阵的基本运算
矩阵的加减法和数乘
矩阵的加减法和数乘其实都非常好理解。
对于两个行数和列数分别相等的矩阵 \(A\) 和 \(B\)，定义它们的和为 \(A+B\)，\(A+B\) 也是一个矩阵，大小和 \(A\) 与 \(B\) 也相等，其元素为原来两个矩阵的对应元素相加。
若：
那么
另外，对于矩阵 \(A\) 和实数 \(k\in \mathbb R\)，定义矩阵 \(kA\) 为矩阵 \(A\) 的数乘，其中每一个元素均为原矩阵对应元素的 \(k\) 倍。
若：
则：
矩阵乘法
对于一个 \(m\times n\) 的矩阵 \(A\)，和一个 \(n\times p\) 的矩阵 \(B\)，我们定义它们可以相乘，相乘的结果为 \(A\times B\)，为一个 \(m\times p\) 的矩阵。
假如 \(A\) 的第 \(i\) 行 \(j\) 列的元素为 \(a_{ij}\)，\(B\) 的第 \(i\) 行 \(j\) 列 的元素为 \(b_{ij}\)，\(A\times B\) 的第 \(i\) 行 \(j\) 列的元素为 \(c_{ij}\)，那么有：
下面的这个矩阵乘法或许可以作为一个例子：
当然，从上面这段描述里，我们也可以证明矩阵乘法的一些性质：
- 两个矩阵的乘法不满足交换律。对矩阵 \(A\) 和 \(B\)，\(A\times B\) 有定义，\(B\times A\) 不一定有定义。
- 两个矩阵可以相乘，当且仅当前一个矩阵的列数等于后一个矩阵的行数。
- 矩阵的乘法满足结合律，即 \(A\times B\times C=A\times (B\times C)\)。
那么，就像 \(1\) 与任何数相乘都得那个数本身一样，是否存在一个矩阵 \(I\)，使得任意的一个 \(m\) 行 \(n\) 列的矩阵与 \(I\) 相乘，都得那个矩阵本身呢？
或者，用形式化的语言讲述：求满足条件的矩阵 \(I\in \mathbb{R}^{n\times n}\)，使得 \(\forall A\in \mathbb{R}^{m\times n}\)，都有 \(A\times I=A\)，其中 \(\mathbb{R}^{m\times n}\) 为所有 \(m\) 行 \(n\) 列的矩阵的集合。
这样的 \(I\) 是存在的！考虑下面这个矩阵：
这个矩阵的对角线上的元素为 \(1\)，其余元素均为 \(0\)，那么对于任意的 \(A\in\mathbb{R}^{m\times n}\)，一定有 \(A\times I=A\)。证明略去。
像这样的，对角线上的元素为为 \(1\)，其余均为 \(0\) 的矩阵，我们称为单位矩阵，它充当矩阵乘法中的单位元。
所以，一个方程组，也可以表示为矩阵乘积的形式，例如：
设
那么原方程组就等价于 \(AX=B\)，两边同时乘以 \(A^{-1}\)，即得 \(X=A^{-1}B\)。
这种办法也可以用来求方程组的解，但是它涉及到矩阵求逆相关的知识，这里暂且不加赘述。
矩阵的转置
将一个矩阵把行和列互换的操作，称为矩阵的转置，记作 \(A^\text T\)，其中 \(A\) 代表要转置的矩阵。
设：
那么：
例如：
行列式
如果一个矩阵的行数和列数相等，那么我们称这个矩阵为一个方阵。
对于一个方阵 \(A\in\mathbb R^{n\times n}\)，我们可以定义它的行列式 \(|A|\)。
设 \(A\) 的第 \(i\) 行 \(j\) 列的元素为 \(a_{ij}\)，一个 \(1\sim n\) 的全排列为 \(c_1, c_2, \cdots, c_n\)，它的逆序对的数目为 \(\tau(c_1c_2\cdots c_n)\)，那么我们定义：
例如，对于一个二阶行列式，有 \(\left|\begin{matrix}a & b\\c & d\end{matrix}\right|=ad-bc\)。
行列式可以通过 Cramer 法则来求解多元一次方程组，我们会在接下来的文章中进行介绍。
Gauss 消元法
到此为止，我们终于介绍完了所有的前置知识，终于可以开始讲 Gauss 消元法了！
要注意的是，Gauss 消元法并非是什么高端的操作，只是一个求解多元方程组的一个可行性方法。
它在实际生活中的应用并不比加减消元广泛，但是由于思维简单，代码清晰而广泛用于计算机求解多元一次方程组的领域。
Gauss 消元法德步骤
它的步骤如下：
- 每次选择一个当前元 \(x_i\)，找到使当前元的系数的绝对值最大的一个方程（在它的增广矩阵矩阵中，这个方程对应着同一行），通过初等行变换，把这一行和第 \(i\) 行交换。
- 把第 \(i\) 行中的当前元的系数，通过初等行变换变为 \(1\)。
- 通过第 \(i\) 行的系数为 \(1\) 的当前元，把所有其他行的当前元的系数都变为 \(0\)。
- 不断重复这个操作，直到每一行都至多有一个 \(1\)，这时，我们就可以得到这个方程组的解，事实上，也有可能无解，或者有无数组解。
『光这么说太抽象了吧？』
那么，我们就以下面这个方程为例，详细地讲一下：
首先，我们可以把它的增广矩阵写出来：
我们发现， \(x_1\) 的系数的绝对值最大的行是第三行，把它与第一行交换，变成了：
把第一行中 \(x_1\) 的系数同时除以 \(3\)，变成 \(1\)，变成：
然后把第一行乘以 \(-1\)，加到第二行，就变成了：
现在，\(x_2\) 系数的绝对值的最大值就在第二行，不用交换了，直接把它变成 \(1\)，变成下面的这个矩阵：
接下来，利用第二行的这个 \(1\)，把另外两行的 \(x_2\) 项全都消去为 \(0\)，变成了：
把第三行乘以一个 \(\dfrac{3}{7}\)，使得 \(x_3\) 项的系数变为 \(1\)，就是：
然后利用这个 \(x_3\)，把另外两行的 \(x_3\) 项消去，就成为了
『好漂亮的矩阵！如果把它变回方程的形式呢？』
那就是：\(\begin{cases}x_1=1,\\x_2=2,\\x_3=3.\end{cases}\)
所以，我们就能通过这种 Gauss 消元法得到方程组的解了。
然而我们会发现，这种做法效率并不高，所以我们平时动笔计算的时候，一般都是用加减消元而非 Gauss 消元。
但是，作为一个具有固定的步骤流程的算法，Gauss 消元比加减消元那种灵活的方法更适合计算机进行了。
代码可以这么写：
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;
const int MAXN = 501;
const double EPS = 1e-5;
int n;
double A[MAXN][MAXN];
int main(void) {
scanf("%d", &n);
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= n; ++j)
scanf("%lf", &A[i][j]);
scanf("%lf", &A[i][n + 1]);
}
//Gauss 消元主体
for(int i = 1; i <= n; ++i) {
int tmp = i; //tmp 表示主元系数绝对值最大的行
for(int j = i + 1; j <= n; ++j)
if(fabs(A[tmp][i]) < fabs(A[j][i])) //打擂台找到绝对值最大的系数
tmp = j;
for(int j = 1; j <= n + 1; ++j)
swap(A[i][j], A[tmp][j]); //将该行与当前行交换
if(fabs(A[i][i]) <= EPS) {
//如果主元的绝对值最大也只是 0，则代表该元
//可以被前面的几个方程消去，矩阵不满秩，则无解
//注意 double 类型不能直接进行比较，可以用 EPS 比较法
printf("No Solution");
return 0;
}
double cur = A[i][i];
for(int j = 1; j <= n + 1; ++j) A[i][j] /= cur; //将该行主元系数化为 1
//枚举每一行，消去该行的主元
for(int j = 1; j <= n; ++j) {
if(j == i) continue;
double div = A[j][i];
for(int k = 1; k <= n + 1; ++k)
A[j][k] -= div * A[i][k];
}
}
for(int i = 1; i <= n; ++i)
printf("%.2f\n", A[i][n + 1]);
return 0;
}
//by CaO
我们可以看出，这个算法的时间复杂度为 \(\mathcal O(n^3)\)。
线性方程组的解的结构定理
正如上面代码中有这么几行：
if(fabs(A[i][i]) <= EPS) {
printf("No Solution");
return 0;
}
我们可以发现，一个方程组不一定是有唯一解的。
考虑下面这个方程组：
这个方程组当然是没有解的，我们如果把它写成增广矩阵的形式，那就是：
我们进行消元之后，可以得到这样的一个矩阵：
我们会发现，最下面这一行，其实代表了一个不可能成立的等式：\(0=-3\)。
如果，再多试几组？
如果手玩一下上面的两个例子，我们会发现，前一个方程是无解的，而后面的这个方程则有无数组解。
第一个方程在进行 Gauss 消元后，得到的矩阵为：
最下面一行代表的方程为 \(0=1\)，显然不可能成立。
而第二个则可以化为：
它的最下面一行呢？代表了一个没有任何意义的等式：\(0=0\)。像这样的行，我们称之为零行。
把它推广到一般形式：对于一个增广矩阵，把它经过初等行变换之后，可以变成像上面的矩阵，我们称之为行阶梯型矩阵。
假定一个有 \(n\) 个未知数方程组，它对应着的增广矩阵，消元后的行阶梯型矩阵中，非零行的数量不大于 \(n\)，那么这个方程组就没有唯一解。
另外，如果这个方程的增广矩阵中，有一个类似这样的行：
其中 \(k\neq0\) 的话，这个方程就是无解的。
我们定义这个矩阵经过消元后的行阶梯型矩阵矩阵中，非零行的数量为这个矩阵的秩，记作 \(\operatorname{rank}(A)\)，其中 \(A\) 为原来的增广矩阵。
上面的这个命题，也可以这么表述：假如一个方程对应的系数矩阵为 \(A\)，增广矩阵为 \(B\)，那么：
- 方程组有唯一解 \(\Leftrightarrow \operatorname{rank}(A)=\operatorname{rank}(B)=n\)；
- 方程组有无数组解 \(\Leftrightarrow \operatorname{rank}(A)=\operatorname{rank}(B)<n\)；
- 方程组无解 \(\Leftrightarrow \operatorname{rank}(A)\neq\operatorname{rank}(B)\)。
Gauss 消元法的其他应用
矩阵求逆
前面说过，对于一个线性方程组，我们可以把它写成矩阵乘法的形式。例如下面这个线性方程组：
我们可以把它写成 \(AX=B\) 的形式，其中：
所以，把 \(AX=B\) 做一个变形，就得到了 \(X=A^{-1}B\)，其中 \(A^{-1}A=I\)，\(I\) 是 \(n\) 阶单位矩阵。
『所以，如果把 \(A^{-1}\) 求出来，就可以得到所有的解了！那么，怎么求这个 \(A^{-1}\) 呢？』
好吧，求矩阵的逆元不是一件简单的事，但是我们对此并不是束手无策的。
我们依旧可以通过 Gauss 消元法，求出 \(A\) 的逆矩阵。
求逆矩阵的思路
下面关于矩阵求逆的问题，都以方阵为主题展开，即我们暂且不谈行数与列数不等的矩阵。
首先，我们要清楚一点：什么样的矩阵才有逆矩阵？或者说，什么样的矩阵才是可逆的？
前人通过证明，发现，对于一个 \(n\times n\) 的方阵 \(A\)，当且仅当 \(\operatorname{rank}(A)=n\) 时，这个矩阵可逆。
既然上文，我们定义了 \(\operatorname{rank}\) 指的是 Gauss 消元后，这个矩阵的非零行的数目，那么一个很自然的想法就是把它通过 Gauss 消元，简化为行阶梯型矩阵。
现在给定一个两个方阵：
我们定义左乘操作 \(A\mid B\) 为把它们左右拼接在一起，即：
那么对于方阵 \(A\mid I\)，其中 \(I\) 是单位矩阵。可以证明，经过 Gauss 消元法，把左侧变为单位矩阵之后，右侧的 \(I\) 也会相应地变成原来 \(A\) 的逆矩阵 \(A^{-1}\)。即：\(A\mid I \stackrel{\text{Gauss 消元}}\rightarrow I\mid A^{-1}\)。
如果 \(A\) 不能消元成为单位矩阵，那就意味着 \(A\) 不是满秩的，所以就不存在逆矩阵。
这也就是说，通过 Gauss 消元法，我们可以得到一个矩阵的逆矩阵了，但它的时间复杂度依旧是 \(\mathcal{O}(n^3)\)。
代码
代码如下，因为这个代码是针对洛谷上的模板题 P4783 【模板】矩阵求逆给出的，所以用了对 \(10^9+7\) 取模而非使用 double：
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
const int MAXN = 1001;
const ll MOD = 1000000007;
ll inv(ll x) { //Fermat 小定理求逆元
ll ans = 1, p = MOD - 2;
while(p) {
if(p & 1) ans = (ans * x) % MOD;
x = (x * x) % MOD; p >>= 1;
}
return ans;
}
int n;
ll A[MAXN][MAXN];
int main(void) {
scanf("%d", &n);
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= n; ++j)
scanf("%lld", &A[i][j]);
A[i][n + i] = 1;
}
//Gauss 消元求矩阵的逆
for(int i = 1; i <= n; ++i) {
int tmp = i;
for(int j = i + 1; j <= (n << 1); ++j)
if(A[tmp][i] < A[j][i]) //打擂台找到绝对值最大的系数
tmp = j;
//判定矩阵不可逆，不解释
if(A[i][i] == 0) {
printf("No Solution");
return 0;
}
ll cur = inv(A[i][i]);
for(int j = 1; j <= (n << 1); ++j)
A[i][j] = (A[i][j] * cur) % MOD; //该行主元系数化为 1
//和 Gauss 消元一毛一样的套路 =v=
for(int j = 1; j <= n; ++j) {
if(i == j) continue;
ll div = A[j][i];
for(int k = 1; k <= (n << 1); ++k)
A[j][k] = ((A[j][k] - div * A[i][k]) % MOD + MOD) % MOD;
}
}
for(int i = 1; i <= n; ++i) {
for(int j = 1; j <= n; ++j)
printf("%lld ", A[i][n + j]);
printf("\n");
}
return 0;
}
//by CaO
时间复杂度分析
既然矩阵求逆的算法是基于 Gauss 消元法的，那么它的时间复杂度自然与 Gauss 消元法相同，都是 \(\mathcal O(n^3)\)。
行列式求值
除了 Gauss 消元之外，我们还有另外一种方法来求 \(n\) 元线性方程组的解— —Cramer 法则。
Cramer 法则：给定一个 \(n\) 元线性方程组，其未知数分别记作 \(x_1,x_2\cdots,x_n\)。
设其矩阵表示形式为 \(AX=B\)，其中 \(A\) 为该方程组的增广矩阵，\(X\) 为由未知数组成的向量，\(B\) 为由常数项组成的向量，\(A_i\) 为 \(A\) 中的第 \(i\) 列被 \(B\) 取代后的矩阵，若 \(A\ne 0\)，则 \(x_i=\dfrac{|A_i|}{|A|}\)；否则该方程组无解。
『这样说，如果能够求出这些矩阵对应行列式，或许就可以求出方程组的解了？』
说的没错，所以，我们应该如何求一个行列式的值呢？
对于一个 \(n\) 阶行列式 \(|A|\)，我们当然可以用定义来求解它的值：
暴力枚举每个 \(1\sim n\) 的全排列，用 cdq 分治求逆序对，复杂度是 \(\Theta(n!\cdot n\log n)\)。这样的复杂度当然是要爆炸的了 QwQ。
所以我们还是应该研究研究行列式的一些性质吧。
经过前人的研究，我们发现了行列式的一下这些性质：
- 假如一个 \(n\) 阶矩阵 \(A\) 是上三角矩阵，即对于任意的 \(i>j\)，\(a_{ij}=0\)，那么 \(|A|=\prod_{i=1}^{n}a_{ii}\)；
- 对于矩阵 \(A\)，将其中任意一行的每个元素同时乘以一个常数 \(k\)，得到矩阵 \(A'\)，则 \(|A'|=k|A|\)；
- 一个矩阵的行列式等于其转置的行列式，即 \(|A|=|A^\text T|\)；
- 对于矩阵 \(A\)，将其中任意两行互换之后，得到矩阵 \(A'\)，则 \(|A|=-|A'|\)；
- 对于矩阵 \(A\)，将其中一行加到另一行上，得到矩阵 \(A'\)，则 \(|A|=|A'|\)。
以上五条性质，证明留给读者作为练习。
我们可以发现，性质 \(2,4,5\) 分别对应了矩阵的三个初等行变换的过程。所以，基于初等行变换的 Gauss 消元法，可以把矩阵简化为一个上三角矩阵，然后求值。
下面是洛谷 P7112【模板】行列式求值的代码，在这道题的题目要求中，要求把行列式的值对一个数 \(p\) 取模，如果 \(p\) 不是个质数的话，就要通过更相减损来消元，而不能使用逆元了。
#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int MAXN = 1001;
int n;
ll p, ans = 1, A[MAXN][MAXN];
int main(void) {
scanf("%d%lld", &n, &p);
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) {
for(int j = i + 1; j <= n; ++j) {
while(A[i][i]) {
ll div = A[j][i] / A[i][i]; //辗转相除约去第 j 行第 i 列的数
for(int k = 1; k <= n; ++k)
A[j][k] = ((A[j][k] - div * A[i][k]) % p + p) % p;
for(int k = 1; k <= n; ++k)
swap(A[i][k], A[j][k]);
ans *= -1; //注意换行的时候，行列式要反号
}
for(int k = 1; k <= n; ++k) //把第 j 行还原回去
swap(A[i][k], A[j][k]);
ans *= -1;
}
}
for(int i = 1; i <= n; ++i)
ans = ((ans * A[i][i]) % p + p) % p; //找到对角线
printf("%lld", ans);
return 0;
}
//by CaO
当然，可以证明，均摊的时间复杂度为 \(\mathcal O(n^2(n+\log p))\)。
例题
本题目列表会持续更新。
