高斯消元学习笔记

· · 算法·理论

对于如下一方程组：

\begin{cases} 2x+3y=7\\ 4x-2y=6 \end{cases}

可以使用加减消元法，易得出解：

\begin{cases} x=2\\ y=1 \end{cases}

我们将方程组扩展一下，变成一个 n 元一次方程组：

\begin{cases} A_{1,1}x_1+A_{1,2}x_2+A_{1,3}x_3+\cdots+A_{1,n}x_n=B_1\\ A_{2,1}x_1+A_{2,2}x_2+A_{2,3}x_3+\cdots+A_{2,n}x_n=B_2\\ A_{3,1}x_1+A_{3,2}x_2+A_{3,3}x_3+\cdots+A_{3,n}x_n=B_3\\ \cdots\\ A_{n,1}x_1+A_{n,2}x_2+A_{n,3}x_3+\cdots+A_{n,n}x_n=B_n \end{cases}

称如上方程组为线性方程组，采用同样的消元方法，称为高斯消元。

高斯消元

以开篇的方程组为例。

我们将所有的系数提出来，得到如下系数矩阵：

\begin{pmatrix} 2 & 3\\ 4 & -2 \end{pmatrix}

将右侧的常数提出来，得到常数矩阵：

\begin{pmatrix} 7\\ 6 \end{pmatrix}

合并一下，就得到了增广矩阵：

\left(\begin{array}{cc|c} 2 & 3 & 7 \\ 4 & -2 & 6 \end{array}\right)

这就是我们所要操作的矩阵。

幻想一下，如果现在我们能得到：

\begin{cases} a_1x+0y=b_1\\ 0x+a_2y=b_2 \end{cases}

的方程，我们是不是就能直接求出解，对吧：

\begin{cases} 1x+0y=\dfrac{b_1}{a_1}\\ 0x+1y=\dfrac{b_2}{a_2} \end{cases}

转换为矩阵的话，就是如下的样子：

\begin{pmatrix} 1 & 0 & \dfrac{b_1}{a_1}\\ 0 & 1 & \dfrac{b_2}{a_2} \end{pmatrix}

替换一下：

\begin{pmatrix} 1 & 0 & x\\ 0 & 1 & y \end{pmatrix}

这就是我们最终所期望得到的矩阵。如果元数变大，那么整体的 1 也会呈现一个阶梯形的样子，所以我们称这种矩阵为简化阶梯形矩阵。

特别的，若最终矩阵存在系数全为 0 但常数不为 0 的一行，则说明这个方程组无解。若存在一行数字全为 0，那么则说明这个方程组有无数解。特别的，判断无解的优先级大于判断无数解。

比如，最终矩阵如下：

\begin{pmatrix} 1 & 0 & 0 & 4\\ 0 & 1 & 0 & 3\\ 0 & 0 & 0 & 0 \end{pmatrix}

这意味着对于任意值的 x_3，原方程组都成立，所以我们称 x_3 为自由元，x_1，x_2 为主元。

综上，我们要我们的增广矩阵向简化阶梯形矩阵靠拢，接下来就要用到高斯消元所需要的三类操作了。

如果是一个正常方程组，使用加减消元法总需要拿一个方程的若干倍去减另一个方程，那么这在高斯消元里对应的操作就是：将其中一行的若干倍（不为 0 ）加到另一行上。

在中间，我们会有一个将一个方程变为原来若干倍的操作，对应在高斯消元里就是：一行所有数乘以一个非零数字。

还有一个操作就是：交换某两行的数字。

我们就需要这三个操作来得到最终的矩阵，当然这三个操作也成初等行变换。

这里我们就不举例了，读者可以自行下去推推。

接下来考虑程序实现。

首先，我们每组操作的目的是确定一个元，我们第 i 个元用第 i 行。

先看看如下方程组：

\begin{cases} 0.000001x+y=1\\ x+y=2 \end{cases}

如果我们选择方程 1 当主元，那么在实际程序中我们可能会丢失大量精度，不如用方程 2 当主元。

所以我们在每组操作里，首先要做的是：保证 A_{i,i} 不为 0，并找到 A_{j,i} 不为 0 且绝对值最大的第 j 列，并交换两行。

接下来，我们就要开始消元了。

对于第 i 列，我们要让除 A_{i,i} 以外都变为 0，具体的来说，就是对于所有的 j\ (j\ne i) ，第 j 行减去 \dfrac{A_{j,i}}{A_{i,i}}。

当 1\sim N 行都处理完了过后，我们需要让所有的 A_{i,i} 变为 1，即 对于第 i 行，整行除以 A_{i,i}，得到简化阶梯形矩阵，最后一列就是方程的解。

具体代码实现如下：

:::success[代码实现]

时间复杂度为 O(n^3)。

const double EPS=1e-9;//精度是个好问题
int gauss() {
    int rank = 0;

    for (int col = 1; col <= n; col++) {
        int ma = rank + 1;
        for (int i = rank + 1; i <= n; i++) {
            if (fabs(a[i][col]) > fabs(a[ma][col]))
                ma = i;
        }

        if (fabs(a[ma][col]) < EPS) continue;

        for (int j = col; j <= n + 1; j++)
            swap(a[rank + 1][j], a[ma][j]);

        for (int i = 1; i <= n; i++) {
            if (i == rank + 1) continue;
            double rate = a[i][col] / a[rank + 1][col];
            if (fabs(rate) < EPS) continue;
            for (int j = col; j <= n + 1; j++)
                a[i][j] -= a[rank + 1][j] * rate;
        }

        rank++;
    }

    for (int i = rank + 1; i <= n; i++) {
        bool all_zero = true;
        for (int j = 1; j <= n; j++) {
            if (fabs(a[i][j]) > EPS) {
                all_zero = false;
                break;
            }
        }
        if (all_zero && fabs(a[i][n + 1]) > EPS)
            return -1;//无解
    }

    if (rank < n) return 1;//无数解

    for (int i = 1; i <= n; i++) {
        if (fabs(a[i][i]) < EPS) continue;
        for (int j = n + 1; j >= i; j--)
            a[i][j] /= a[i][i];
    }

    return 0;//唯一解
}

:::

习题

:::info[「SDOI2006」线性方程组] 模板题。借用此题熟悉高斯消元。 :::

异或方程组

如果 \forall 1\le i\le N,\ x_i\in\{0,1\}，那么存在异或方程组如下：

\begin{cases} A_{1,1}x_1\oplus A_{1,2}x_2\oplus A_{1,3}x_3\oplus \cdots\oplus A_{1,n}x_n=B_1\\ A_{2,1}x_1\oplus A_{2,2}x_2\oplus A_{2,3}x_3\oplus \cdots\oplus A_{2,n}x_n=B_2\\ A_{3,1}x_1\oplus A_{3,2}x_2\oplus A_{3,3}x_3\oplus \cdots\oplus A_{3,n}x_n=B_3\\ \cdots\\ A_{n,1}x_1\oplus A_{n,2}x_2\oplus A_{n,3}x_3\oplus \cdots\oplus A_{n,n}x_n=B_n \end{cases}

我们需要求解，同样可以采用高斯消元。

已知，异或运算可以看做不进位的加法。我们仍然可以写出增广矩阵，在执行高斯消元的过程中，把加减换为异或，且不用执行乘除法。

同上，我们可以写出过程：

若最后存在一行系数全为 0 但是常数为 1，那么方程组无解；如果有一行全为 0，那么该行的 x_i 为自由元。有解情况下，如果自由元一共有 c 个，那么最后解的数量为 2^c 个。

在 n 小的时候，我们可以使用 int 和 long long 来存储，否则就可以使用 bitset。读者可以自行下去学习 bitset。

习题

:::info[开关问题] 注意到每个开关只有开（1）或关（0）两种状态，所以我们可以设第 i 个开关的状态为 x_i。列出 N 个方程。

对于系数，如果一个开关 j 会影响开关 i 则那么，令 A_{i,j}=1，表示会影响。因为情况较少，读者可以自己下去枚举验证。

令第 i 个开关开始时为 s_i，结束时为 e_i。则那么以得出最终第 i 行的常数为 s_i\oplus e_i。

这样我们易列出：

\begin{cases} A_{1,1}x_1\oplus A_{1,2}x_2\oplus A_{1,3}x_3\oplus \cdots\oplus A_{1,n}x_n=s_1\oplus e_1\\ A_{2,1}x_1\oplus A_{2,2}x_2\oplus A_{2,3}x_3\oplus \cdots\oplus A_{2,n}x_n=s_2\oplus e_2\\ A_{3,1}x_1\oplus A_{3,2}x_2\oplus A_{3,3}x_3\oplus \cdots\oplus A_{3,n}x_n=s_3\oplus e_3\\ \cdots\\ A_{n,1}x_1\oplus A_{n,2}x_2\oplus A_{n,3}x_3\oplus \cdots\oplus A_{n,n}x_n=s_n\oplus e_n \end{cases}

高斯消元即可。 :::

线性空间

因为我们老师说这和高斯消元有关所以也讲了，给一位初一生干懵了。

:::warning[如果你不懂向量的话请看这里] 向量（Vector），顾名思义就是有方向的量，一个 k 维向量对于坐标系可以用 (a_1,a_2,\dots,a_k) 表示。画在图上的话，就是一个带箭头的线段，线段的方向表示向量的方向，线段的长度表示向量的大小。

与向量对应的就是数量（标量），标量只有大小，没有方向。我们这里所讨论的标量都是实数。

向量也有运算，但在这里我们只用向量加法（a+b，其中 a,b 是向量）与标量乘法（k\times a，其中 k 是实数标量，a 是向量）。 :::

在一个非空向量集合 V 中,其上定义了向量加法和标量乘法，可以通过集合中的已有的向量进行运算后所得到的向量仍在集合 V 中，称这样一个系统为线性空间。

对于一组向量 \alpha_1,\alpha_2,\dots,\alpha_n，若存在不全为零的标量 k_1,k_2,\dots,k_n，使得：

k_1\alpha_1+k_2\alpha_2+\cdots+k_n\alpha_n=0

则称这组向量线性相关；否则称线性无关。

若线性空间 V 中的一组向量 \alpha_1,\alpha_2,\dots,\alpha_n 满足：

  1. 线性无关；

则称这组向量为 V 的基（基底）。

举个例子。

考虑二维平面上的所有向量组成的线性空间 \mathbb{R}^2。

发现所有的向量 (x,y) 都能由 \{(1,0),\ (0,1)\} 组合出（x(1,0)+y(0,1)）。所以 \{(1,0),\ (0,1)\} 是该线性空间的基。

其实，我们可以把矩阵的每一行看成一个 M 维向量，矩阵的 n 个向量通过向量运算能表示出来的所有向量构成一个线性空间。我们将这个矩阵进行高斯消元，得到的简化阶梯形矩阵的所有非 0 行向量线性无关。我们执行的三种操作就是标量乘法（只会乘以非 0 向量）和向量加法，并不会改变这个线性空间。所以，我们所得到的所有非 0 行向量就是这个线性空间的基，个数称为这个线性空间的秩。

习题

::::info[[JLOI2015] 装备购买]

转换一下题意，易得出我们所需要的向量应该是线性无关的，所以建立矩阵，进行高斯消元。

对于要购买尽可能多的装备，我们可以使用贪心。在选择行时，优先选择费用较少的。

贪心正确性证明如下：

:::info[证明]{open}

把所有的 z 按价格从小到大排序。

若 z_k 与 z_1, \dots, z_{k-1} 中已选的向量线性无关，则存在一个最优基包含 z_k。

证明：假设最优基 G 不含 z_k。因 G 是基，z_k 可由 G 表示。又因 z_k 不能被前面已选向量表示，所以表示式中必用到 G 中某个“非前面已选”的向量 z_p。移项可得 z_p 能被 \{z_k\} \cup (G - \{z_p\}) 表示，故：

G' = \{z_k\} \cup (G - \{z_p\})

与 G 张成同一空间，G' 也是基。而 z_k 比 z_p 便宜（k < p），所以 G' 花费更低，矛盾。

故贪心每次选的向量都能扩展成一个最优基，最终得到的就是最优基。贪心正确。

::: ::::

后记

在写完全文初稿之后，使用过 AI 查错，但保证个人贡献大于 AI 贡献。

全文 6.9k+ 字，完。