高斯消元学习笔记

· · 算法·理论

距离 NOIP 还剩两个多星期，写篇文章攒攒 RP，同时巩固一下高斯消元。

引入

高斯消元是求解线性方程组的经典算法，由德国数学家高斯提出。

例一（P2455 线性方程组）

这道题目就是高斯消元的模板。

题意

解以下 n (1 \le n \le 50) 元线性一次方程组：

\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}

其中 a_{1, 1}, \cdots ,a_{n, n}，b_1, \cdots, b_n 均为给定常数。

若存在唯一解，则输出；若无解，输出 -1；若存在无穷多实数解，输出 0。

高斯消元

消元法

首先回顾一下初中学过的解二/三元一次方程组的方法。主要思想便是消元。而消元方法又分为两种：代入消元法和加减消元法。

前者顾名思义，将一未知数用含有其它未知数的代数式表示，并将其带入到别的方程中，从而消去这一未知数。比如对于下面的方程组：

\begin{cases} 4x + y = 100 \\ x - y = 100 \end{cases}

由第二个方程可得

x = 100 + y

将其带入方程组中第一个方程，消元 x 得

400 + 5y = 100

解得

y = -60

再将其回代到方程 x = 100 + y 得

x = 40

再说后者，即将方程组中的一方程倍乘某个常数加到另外的方程中去，实现消元。仍以上面的方程组为例：

\begin{cases} 4x + y = 100 \\ x - y = 100 \end{cases}

直接把两方程相加，消元 y 得

5x = 200

解得

x = 40

再把 x = 40 带入到方程组中第二个方程可得

y = -60

通过比较我们发现：如果不考虑回代，那么显然加减消元法更方便实现（因为只要令每一项系数相加），所以我们今天重点考虑加减消元法。

我们总结出：对方程组做初等变换，解不变。具体来说，我们可以：

高斯消元法

首先，考虑到消元过程中参与计算和发生改变的是方程中各变量的系数，和各变量无关，方便起见，我们把题目中的线性方程组的系数表示成一个 n \times n 矩阵：

A = \begin{bmatrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{bmatrix}

对于方程右边的常数，作为一个 n \times 1 列向量 b ，放在矩阵 A 的右边，得到 A 的增广矩阵：

[A | b] = \left [ \begin{array} {c:c} \begin{matrix} a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\ \end{matrix} & \begin{matrix} b_1 \\ b_2 \\ \vdots \\ b_n \\ \end{matrix} \end{array} \right ]

此时线性方程组的初等变换对应矩阵的初等行变换（列可不能乱动）。

通过初等行变换，把增广矩阵化为行阶梯矩阵，模拟消元过程，从而求解线性方程组，这就是最基础的高斯消元法。

到这里，你可能有些没听懂（好吧其实说的是笔者自己），所以我们举个简单的例子：

解以下方程组：

\begin{cases} x + 2y + z = 7 \\ 2x - y + 3z = 7 \\ 3x + y + 2z = 18 \\ \end{cases}

先写出该方程组对应的的增广矩阵：

\left [ \begin{array} {c:c} \begin{matrix} 1 & 2 & 1 \\ 2 & -1 & 3 \\ 3 & 1 & 2 \\ \end{matrix} & \begin{matrix} 7 \\ 7 \\ 18 \\ \end{matrix} \end{array} \right ]

接下来化为行阶梯矩阵：

\left [ \begin{array} {c:c} \begin{matrix} 1 & 2 & 1 \\ 2 & -1 & 3 \\ 3 & 1 & 2 \\ \end{matrix} & \begin{matrix} 7 \\ 7 \\ 18 \\ \end{matrix} \end{array} \right ] \xrightarrow{r_2 - 2r_1} \left [ \begin{array} {c:c} \begin{matrix} 1 & 2 & 1 \\ 0 & -5 & 1 \\ 3 & 1 & 2 \\ \end{matrix} & \begin{matrix} 7 \\ -7 \\ 18 \\ \end{matrix} \end{array} \right ] \\ \xrightarrow{r_3 - 3r_1} \left [ \begin{array} {c:c} \begin{matrix} 1 & 2 & 1 \\ 0 & -5 & 1 \\ 0 & -5 & -1 \\ \end{matrix} & \begin{matrix} 7 \\ -7 \\ -3 \\ \end{matrix} \end{array} \right ] \xrightarrow{r_3 - r_2} \left [ \begin{array} {c:c} \begin{matrix} 1 & 2 & 1 \\ 0 & -5 & 1 \\ 0 & 0 & -2 \\ \end{matrix} & \begin{matrix} 7 \\ -7 \\ 4 \\ \end{matrix} \end{array} \right ]

发现有 -2z = 4 ，解得 z = -2 。

接下来自底向上依次回代可得 y = 1 ， x = 7 。

通过上面的例子，大致流程可以轻松看出：枚举列 j = 1,2,\cdots,n，用第 j 行将第 j+1 到第 n 行中的第 j 列消去。

具体来说，对于每个 j 枚举 i = j+1,j+2, \cdots ,n ，令 r_i - \frac{a_{i, j}}{a_{j, j}} r_j，即可消去第 j+1 到第 n 个方程中的未知数 x_j 。

理想情况下矩阵变为上三角矩阵。

最后做自底向上的回代即可。

实现

方便起见，定义初等行变换函数，用来实现形如 r_1+kr_2 的变换。

void row_add(double r1[], double r2[], double k) {
  for (int j = 1; j <= n + 1; j++) r1[j] += k * r2[j];
}

下面考虑增广矩阵化为行阶梯矩阵。之前提到，枚举列 j，对于每个 j 枚举 i = j+1,j+2,\cdots,n ，令 r_i - \frac{a_{i, j}}{a_{j, j}} r_j。但是一种特殊情况为 a_{j, j} = 0 ，此时变换无意义。因此我们需要找到一行使得 x_j 的系数不为零，并且用这一行与第 j 行交换。若对于任意行 r (j \le r \le n) 都有 a_{r, j} = 0 ，那么直接跳过这一列，x_j 成为自由元。注意，为了最终将矩阵化为行阶梯矩阵，下次消元时我们仍从这一行开始寻找主元系数不为零的行，这也意味着对于一列 j，我们并不总是恰好用第 j 行消元，而是需要保存一个变量来记录当前准备用来消元的行，只有消元成功才会把该变量更新成下一行。

int cur = 1;
for (int j = 1; j <= n; j++) {
  int k = cur;
  for (int i = cur + 1; i <= n; i++)
    if (abs(a[i][j]) > abs(a[k][j])) k = i;
    if (abs(a[k][j]) < EPS)
    continue;
  swap(a[cur], a[k]);
  for (int i = j + 1; i <= n; i++)
    row_add(a[i], a[cur], -a[i][j] / a[cur][j]);
  cur++;
}

这里不用第一个主元系数非零的行而是用主元系数最大的行，目的是最小化 \frac{a[i][j]}{a[cur][j]} 从而减小 double 带来的误差。另外浮点数判断是否为零也不能直接判断是否等于零，而是需要用一个极小值（这里 EPS = 10 ^ {-6}）和该数的绝对值作比较。

最后考虑回代。假设数据存在唯一解，经过高斯消元后我们得到一个上三角矩阵。矩阵第 i 行所对应的方程为：

a_{i, i}x_i + a_{i, i+1}x_{i+1} + \cdots + a_{i, n}x_n = b_i

解得

x_i = \frac{b_i - a_{i, i+1}x_{i+1} - \cdots - a_{i, n}x_n}{a_{i, i}}

显然 x_i 只依赖于 x_{i+1} 到 x_{n}。矩阵自底向上回代即可。

for (int i = n; i >= 1; i--) {
  x[i] = a[i][n + 1];
  for (int j = i + 1; j <= n; j++)
    x[i] -= a[i][j] * x[j];
  x[i] /= a[i][i];
}

注意，这么做的前提是存在唯一解，如果无解或有无穷多组解则矩阵并不是上三角矩阵，上述代码第 5 行可能会出现除以 0 的情况，导致 \colorbox{purple}{\textcolor{white}{RE}} 。

到这里我们已经可以解决有唯一解的数据。

判断无解、无穷多组解

首先如果没有出现自由元，则一定存在唯一解。

当出现自由元时，矩阵底部必定存在系数部分全为 0 的行。

在这些行中，若有至少一行右侧常数不为 0，则代表有 0=k (k \neq 0) 的式子，显然无解。

否则如果右侧常数也全为 0，则说明本质不同的方程不足 n 条，此时必定存在无穷多组解。

实现
int flag = 1;
if (--cur < n) {
  for (int i = cur + 1; i <= n; i++)
    if (abs(a[i][n + 1]) > EPS) {
      flag = -1; // 无解
      break;
    }
  if (flag == 1) flag = 0; // 无穷多组解
}

注意，cur 为当前准备用来消元的行，实际消元数量为 cur-1。所以若 cur - 1 < n 则代表有自由元；否则 cur - 1 = n 则代表存在唯一解。

至此，你已经学会了基础的高斯消元法。

高斯-约旦消元法

不过遗憾地告诉你，刚才的基础款高斯消元法我们平时根本不用。因为我们有升级款——高斯-约旦消元法。

注意到上述高斯消元法旨在将增广矩阵变换为行阶梯矩阵，最后再回代。因而我们每次只对当前行之后的行进行消元。那么我们为什么不把当前行上面的行也顺带消元呢？

这样一来，理想情况下矩阵的系数部分最终变为对角矩阵，也就是说每行都是一个一元一次方程，两边同除以系数即可得到方程的解，省去了回代的步骤。

当用第 cur 行消元 x_j 时，直接把系数化为 1。理想情况下，矩阵系数部分最终变为单位矩阵，矩阵最右列即方程组的解，一步到位。

至于无解、无穷多组解的情况，矩阵的系数部分只能化为行最简形。

实现

方便起见定义初等行变换函数，实现形如 r \times \frac{1}{k} 的变换。

void row_div(double r[], double k) {
  for (int j = 1; j <= n + 1; j++) r[j] /= k;
}

易得消元部分代码。

int cur = 1;
for (int j = 1; j <= n; j++) {
  int k = cur;
  for (int i = cur + 1; i <= n; i++)
    if (abs(a[i][j]) > abs(a[k][j])) k = i;
  if (abs(a[k][j]) < EPS)
    continue;
  swap(a[cur], a[k]);
  row_div(a[cur], a[cur][j]);
  for (int i = 1; i <= n; i++)
    if (i != cur) row_op(a[i], a[cur], -a[i][j]);
  cur++;
}

无解、无穷多组解的判定与普通高斯消元法相同，此处略。

时间复杂度

显然是未知数数量的平方乘以方程条数，对于这道题是 O(n^3) 。

例二（P2447 外星千足虫）

题意

有 n (n \le 1000) 个未知数 x_1,x_2,\cdots,x_n，每个未知数的取值不是 0 就是 1。

还有 m (m \le 2000) 个方程联立的方程组：

\begin{cases} a_{1, 1} x_1 \oplus a_{1, 2} x_2 \oplus \cdots \oplus a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{2, n} x_n = b_2 \\ \cdots \\ a_{m,1} x_1 \oplus a_{m, 2} x_2 \oplus \cdots \oplus a_{m, n} x_n = b_m \end{cases}

其中 a_{i, j},b_i \in \{ 0, 1 \} 。

如果存在唯一解，则输出最小的 k，满足只需要根据输入的第 1 到第 k 条方程即可确定唯一解，接着输出所求出的唯一解；

如果存在无穷多组解输出 Cannot Determine；

数据保证不会无解。

高斯消元求解异或方程组

我们发现异或满足交换律、结合律，且逆运算仍为异或，所以我们可以通过把两个方程异或在一起来消元，进而将增广矩阵系数部分化为行最简形。

bitset 优化

由于矩阵为 01 矩阵，把两行异或在一起这一操作不值得线性枚举，可以使用 std::bitset 压位，直接异或，从而把时间复杂度优化到 O(\frac{n^2m}{\omega}) ，其中 \omega 一般为 32。

求最少到第几条方程就可以确定唯一解

一个暴力的方式是直接二分。但是这样不仅麻烦，而且时间复杂度会多一层 \log，故这里不再考虑。

有一个简单的思路：记 id_i 表示矩阵第 i 行是输入的第几条方程，对于枚举的主元，每次从当前行往下找 id 最小的且主元系数非零的行来消元。答案即为 \max \limits_{1 \le i \le n} id_i。

不过事实上只需要找第一个主元系数非零的行即可。原因也很简单，设有行 r_1, r_2 主元系数都不为零，分别为输入的第 id(r_1),id(r_2) 条方程，r_1 在 r_2 上方且 id(r_1) > id(r_2)（也就是竖着的逆序对）。由于初始时按照 id 从小到大排列，所以 r_2 初始时一定在 r_1 上方，当前沦落到其下方必定是因为和原来的这一行发生了交换，这意味着我们之前已经把原来的 r_2 这一行用来消元，那么这一行的 id 一定对答案产生了贡献。又因为这一行的 id 一定比 id(r_1) 和 id(r_2) 都大，所以此时选 r_1 还是 r_2 都不会对答案产生影响。

这个方法同样适用于线性方程组。

参考代码

#include <iostream>
#include <bitset>
using namespace std;

const int N = 1000 + 1, M = 2000 + 5;
bitset<N> a[M];
int n, m, id[M];

void gauss() {
  int cur = 1, mx = 0;
  for (int j = 0; j < n; j++) {
    int k = cur;
    while (k <= m && !a[k][j]) k++;
    if (k > m) {
      cout << "Cannot Determine";
      return;
    }
    swap(a[k], a[cur]); swap(id[k], id[cur]);
    for (int i = 1; i <= m; i++)
      if (i != cur && a[i][j]) a[i] ^= a[cur];
    mx = max(mx, id[cur++]);
  }
  cout << mx << endl;
  for (int i = 1; i <= n; i++)
    cout << (a[i][n] ? "?y7M#\n" : "Earth\n");
}

int main() {
  cin >> n >> m;
  for (int i = 1; i <= m; id[i] = i, i++)
    for (int j = 0; j <= n; j++) {
      char c; cin >> c;
      if (c == '1') a[i][j] = 1;
    }
  gauss();
  return 0;
}

例三（P2962 Lights G）

题意

在 n (1 \le n \le 35) 个点 m (1 \le m \le 595) 条边的无向图中，所有点初始状态均为 0，一次操作可以翻转一个点及其相邻节点的状态（0 变成 1，1 变成 0），求最少的操作次数使得所有点状态均为 1。

数据保证有解。

分析

容易发现，每个点的最终状态即为该点与其相邻节点的操作次数之和的奇偶性，且每个点至多操作一次。

设 x_i 表示是否对点 i 进行一次操作，则有方程组：

\begin{cases} a_{1, 1} x_1 \oplus a_{1, 2} x_2 \oplus \cdots \oplus a_{1, n} x_n = 1 \\ a_{2, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{2, n} x_n = 1 \\ \cdots \\ a_{n,1} x_1 \oplus a_{n, 2} x_2 \oplus \cdots \oplus a_{n, n} x_n = 1 \end{cases}

在第 u 条方程中，a_{u, v} 的取值为 u,v 是否相邻（我们假设每个点都和自己相邻）。

高斯消元解方程组即可。

但是题目只保证有解，并未保证解唯一。因此可能出现自由元的情况。

对自由元的处理

顾名思义，自由元取值不受约束，然而非自由元的取值可能受到自由元的限制。

为使最终取值为 1 的未知数尽可能少，我们只能暴力地搜索每个自由元的取值，没有其他好的办法。

实现

高斯消元时，用 r_j 记录 x_j 作为主元的方程在第几行，若 r_j = 0 则 x_j 为自由元。

由于高斯消元后，矩阵系数部分已经化为行最简形，因此非自由元 x_j 的取值最多依赖 x_{j+1} 到 x_n，所以自底向上搜索即可。

每次搜到一个未知数，若为自由元，分别取值 0,1 继续向上搜；若非自由元， 则根据其作为主元的方程推导出其取值，然后再向上搜。

记得加上最优性剪枝，快得飞起。

参考代码

#include <iostream>
#include <bitset>
using namespace std;

const int N = 35 + 1, M = 595 + 5;
int n, m;
bitset<N> a[M];
int r[N], p[N];
int ans = 2e9;

void gauss() {
  int cur = 1;
  for (int j = 1; j <= n; j++) {
    int k = cur;
    while (k <= m && !a[k][j]) k++;
    if (k > m) continue;
    swap(a[k], a[cur]);
    r[j] = cur;
    for (int i = j + 1; i <= m; i++)
      if (a[i][j]) a[i] ^= a[cur];
    cur++;
  }
}

void dfs(int step, int res) {
  if (res >= ans)
    return;
  if (!step) {
    ans = res;
    return;
  }
  if (!r[step]) {
    p[step] = 0;
    dfs(step - 1, res);
    p[step] = 1;
    dfs(step - 1, res + 1);
    return;
  }
  p[step] = a[ r[step] ][0];
  for (int j = 1; j <= n; j++)
    if (j != step && a[ r[step] ][j])
      p[step] ^= p[j];
  dfs(step - 1, res + p[step]);
}

int main() {
  cin >> n >> m;
  for (int i = 1; i <= n; i++)
    a[i][i] = a[i][0] = 1;
  for (int i = 1, u, v; i <= m; i++) {
    cin >> u >> v;
    a[u][v] = a[v][u] = 1;
  }
  gauss();
  dfs(n, 0);
  cout << ans;
}

例四（P10315 原神，启动！）

题意

有 n 个方碑，每个方碑状态为 0 到 m-1 中的一种（保证 m 为素数）。攻击方碑 i 会使 i 自身与给定的一些方碑状态同时加 1 模 m。

给定每个方碑的初始状态 s_i 和最终状态 t_i，求每个方碑需要攻击的次数 x_i（0 ≤ x_i < m），满足所有方碑从初始状态转换到最终状态。

无解输出 niuza。若有多解，输出任意一组。

分析

和上一题相似，容易得到方程组：

\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n \equiv t_1 - s_1 {\pmod m} \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n \equiv t_2 - s_2 {\pmod m} \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n \equiv t_n - s_n {\pmod m} \end{cases}

不同的是，每条方程都是在模意义下的。

高斯消元求解模意义下的方程组

由同余的性质，仍然可以在同余式两边作乘法，或把两方程相加，唯独不能作除法来使系数化为 1。不过我们可以在同余式两边同乘以主元系数在模 m 意义下的逆元来达到相同的效果。

另外，两式相加时可以给系数和常数都取模。

其余部分与普通的线性方程组无异。

虽然逆元需要每次用快速幂或扩展欧几里得算法来求，但时间复杂度仍为 O(n^3)。

参考代码

#include <iostream>
using namespace std;
using ll = long long;

const int N = 100 + 5;
int n, m, s[N], t[N];
ll a[N][N];
int ans[N];

int qpow(ll a, int b) {
  ll ret = 1;
  while (b) {
    if (b & 1) ret = ret * a % m;
    a = a * a % m;
    b >>= 1;
  }
  return ret;
}

void row_mul(ll r[], int k) {
  for (int j = 1; j <= n + 1; j++) r[j] = (r[j] * k % m + m) % m;
}

void row_add(ll r1[], ll r2[], int k) {
  for (int j = 1; j <= n + 1; j++) (r1[j] += k * r2[j] % m + m) %= m;
}

bool gauss() {
  int cur = 1;
  for (int j = 1; j <= n; j++) {
    for (int i = cur; i <= n; i++)
      if (a[i][j]) {
        swap(a[cur], a[i]);
        break;
      }
    if (!a[cur][j]) continue;
    row_mul(a[cur], qpow(a[cur][j], m - 2));
    for (int i = 1; i <= n; i++)
      if (i != cur) row_add(a[i], a[cur], -a[i][j]);
    cur++;
  }
  for (int i = cur; i <= n; i++)
    if (a[i][n + 1]) return 0;
  for (int i = 1; i <= n; i++)
    for (int j = 1; j <= n; j++)
      if (a[i][j]) {
        ans[j] = a[i][n + 1];
        break;
      }
  return 1;
}

int main() {
  cin >> n >> m;
  for (int i = 1; i <= n; i++) {
    int k; cin >> k;
    a[i][i] = 1;
    while (k--) {
      int j; cin >> j;
      a[j][i] = 1;
    }
  }
  for (int i = 1; i <= n; i++)
    cin >> s[i];
  for (int i = 1; i <= n; i++)
    cin >> t[i], a[i][n + 1] = ((t[i] - s[i]) % m + m) % m;
  bool flag = gauss();
  if (!flag) {
    cout << "niuza";
  } else {
    for (int i = 1; i <= n; i++)
      cout << ans[i] << ' ';
  }
  return 0;
}

总结

高斯还是太强了，%%%。