高斯 - 约当消元法

1|0定义

1|1高斯消元方法

给定 $n$ 元一次方程组

{\begin{cases} a_{1, 1} x_{1} + a_{1, 2} x_{2} + \dots + a_{1, n} x_{n} = b_{1} \\ a_{2, 1} x_{1} + a_{2, 2} x_{2} + \dots + a_{2, n} x_{n} = b_{2} \\ \dots \\ a_{n, 1} x_{1} + a_{n, 2} x_{2} + \dots + a_{n, n} x_{n} = b_{n} \end{cases}

请求出方程组的解的情况：

无解；
无穷多解；
唯一解。

对于这样的问题，我们可以使用 高斯消元法 进行求解，当然高斯消元法有一个回代的过程，代码略长，而且精度较低。

所以我们隆重推出 高斯-约旦消元法 ！！！

回顾一下我们是怎么手算的，一般用的都是 加减消元法，普通高斯和高斯-约旦用的都是加减消元。

在此之前，我们需要了解一下矩阵初等变换。

在线性代数中，矩阵初等行变换 是指以下三种变换类型：

交换矩阵的两行；
用一个非零数 $k$ 乘矩阵的某一行所有元素；
把矩阵的某一行所有元素乘以一个数 $k$ 后加到另一行对应的同一列的元素上；

类似地，把以上的 行 改为 列 便得到 矩阵初等列变换 的定义。

矩阵初等行变换与初等列变换合称为 矩阵初等变换。

若矩阵 $A$ 经过有限次的初等行变换变为矩阵 $B$ ，则矩阵 $A$ 与矩阵 $B$ 行等价；若矩阵 $A$ 经过有限次的初等列变换变为矩阵 $B$ ，则矩阵 $A$ 与矩阵 $B$ 列等价；若矩阵 $A$ 经过有限次的初等变换变为矩阵 $B$ ，则矩阵 $A$ 与矩阵 $B$ 等价。

~~当然列的用不着~~

首先有一个由系数构成的 $n \times n$ 的矩阵

[\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{n, 1} & a_{n, 2} & \dots & a_{n, n} \end{matrix}]

然后是一个由常数构成的 $n \times 1$ 的列向量

[\begin{matrix} b_{1} \\ b_{2} \\ ⋮ \\ b_{n} \end{matrix}]

把它们放在一起构成一个 $n \times (n + 1)$ 的增广矩阵：

[\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, n} & ∣ & b_{1} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, n} & ∣ & b_{2} \\ ⋮ & ⋮ & ⋱ & ⋮ & ∣ & ⋮ \\ a_{n, 1} & a_{n, 2} & \dots & a_{n, n} & ∣ & b_{n} \end{matrix}]

我们遍历每一列，对于第 $i$ 列选出最大的 未处理过的行 上的数作为主元，意味着加减消元后除了主元这一行外其他行的第 $i$ 列都为 $0$ 。

选最大的作为主元的原因是避免精度误差。

然后就是加减消元了。

举个例子

{\begin{cases} 2 x - y + z = 1 \\ 4 x + y - z = 5 \\ x + y + z = 0 \end{cases}

增广矩阵就是

[\begin{matrix} 2 & - 1 & 1 & ∣ & 1 \\ 4 & 1 & - 1 & ∣ & 5 \\ 1 & 1 & 1 & ∣ & 0 \end{matrix}]

先是第 $1$ 列，选 $4$ 为主元。

$4$ 在第 $2$ 行，将第 $2$ 行与第 $1$ 行交换。

[\begin{matrix} 4 & 1 & - 1 & ∣ & 5 \\ 2 & - 1 & 1 & ∣ & 1 \\ 1 & 1 & 1 & ∣ & 0 \end{matrix}]

对于第 $2$ 行，第一列上的 $2$ 是 $4$ 的 $\frac{1}{2}$ 。

[\begin{matrix} 4 & 1 & - 1 & ∣ & 5 \\ 2 - 4 \times \frac{1}{2} & - 1 - 1 \times \frac{1}{2} & 1 - (- 1) \times \frac{1}{2} & ∣ & 1 - 5 \times \frac{1}{2} \\ 1 & 1 & 1 & ∣ & 0 \end{matrix}]

化简得

[\begin{matrix} 4 & 1 & - 1 & ∣ & 5 \\ 0 & - \frac{3}{2} & \frac{3}{2} & ∣ & - \frac{3}{2} \\ 1 & 1 & 1 & ∣ & 0 \end{matrix}]

对第 $3$ 行的处理同理

[\begin{matrix} 4 & 1 & - 1 & ∣ & 5 \\ 0 & - \frac{3}{2} & \frac{3}{2} & ∣ & - \frac{3}{2} \\ 0 & \frac{3}{4} & \frac{5}{4} & ∣ & - \frac{5}{4} \end{matrix}]

现在到了第 $2$ 列，未处理过的是 $2, 3$ 行，选最大的 $\frac{3}{4}$ 为主元。

将第 $3$ 行与第 $2$ 行交换

[\begin{matrix} 4 & 1 & - 1 & ∣ & 5 \\ 0 & \frac{3}{4} & \frac{5}{4} & ∣ & - \frac{5}{4} \\ 0 & - \frac{3}{2} & \frac{3}{2} & ∣ & - \frac{3}{2} \end{matrix}]

消元得

[\begin{matrix} 4 & 0 & - \frac{8}{3} & ∣ & \frac{20}{3} \\ 0 & \frac{3}{4} & \frac{5}{4} & ∣ & - \frac{5}{4} \\ 0 & 0 & 4 & ∣ & - 4 \end{matrix}]

第 $3$ 列，未处理过的是第 $3$ 行，选 $4$ 作主元。

[\begin{matrix} 4 & 0 & 0 & ∣ & 4 \\ 0 & \frac{3}{4} & 0 & ∣ & 0 \\ 0 & 0 & 4 & ∣ & - 4 \end{matrix}]

这样就把原来的矩阵通过初等变换，使得系数矩阵只有主对角线位置的元素非零，即一个对角矩阵。

上面那个矩阵的意思是

{\begin{cases} 4 x = 4 \\ \frac{3}{4} y = 0 \\ 4 z = - 4 \end{cases}

所以再把系数除过去就得到

{\begin{cases} x = 1 \\ y = 0 \\ z = - 1 \end{cases}

请自行检验。

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

当然方程还可能是无解或无穷多解。

考虑一个一元一次方程 $a x = b$ 的解的情况：

无解： $a = 0, b \neq 0$ ；
无穷多解： $a = b = 0$ ；
唯一解： $a \neq 0$ 。
所以当发现某一列的主元为 $0$ 时，因为主元是最大的，所以意味着这一列全都是 $0$ ，那么要么无解，要么有无穷多解。

P3389 【模板】高斯消元法

#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
double eps = 1e-7;
signed main()
{
    int n; cin >> n;
    vector<vector<double>> a(n + 5, vector<double> (n + 6));
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n + 1; j++){
            cin >> a[i][j];
        }
    }
    for(int j = 1; j <= n; j++){ // 枚举列
        int max = j;
        for(int k = j + 1; k <= n; k++){
            if(fabs(a[k][j]) > fabs(a[max][j])){
                max = k; // 得到 j 列的最大系数所在的行
            }
        }
        // 最大元系数调换到第 j 行
        swap(a[j], a[max]);
        if(fabs(a[j][j]) < eps){
            return cout << "No Solution" << '\n', 0;
        }
        // 主元系数变为 1
        for(int k = n + 1; k >= 1; k--){
            a[j][k] = a[j][k] / a[j][j];
        }
        // 其余行系数变成 0 
        for(int i = 1; i <= n; i++){
            if(i == j) continue;
            double tmp = a[i][j];
            for(int k = 1; k <= n + 1; k++){
                a[i][k] = a[i][k] - tmp * a[j][k];
            }
        }
    }
    for(int i = 1; i <= n; i++) printf("%.2lf\n", a[i][n + 1]);
    return 0;
}

1|2关于判断无解和无穷多解

关于判断无解和无穷多解

P2455 [SDOI2006]线性方程组

这下要具体到到底是无解还是无穷多解了。

这就是高斯-约旦的一个缺点：相比于普通高斯，它更难判断无解和无穷多解。

但是也是可以处理的。

具体地，我们用 $r$ 来记录当前行，如果主元为 $0$ ，那么 $continue$ 到下一列，但 $r$ 不变；否则消元后令 $r \leftarrow r + 1$ 。

遍历完所有列后：

若 $r = n$ ，说明有唯一解；
若 $r < n$ ，说明第 $r + 1 \sim n$ 行系数矩阵全都是 $0$ ，若列向量全是 $0$ ，说明有无穷多解，否则无解。

#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
double eps = 1e-7;
signed main()
{
    int n; cin >> n;
    vector<vector<double>> a(n + 5, vector<double> (n + 6));
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n + 1; j++){
            cin >> a[i][j];
        }
    }
    int r = 1;
    for(int j = 1; j <= n; j++){ // 枚举列
        int max = r;
        for(int k = r + 1; k <= n; k++){
            if(fabs(a[k][j]) > fabs(a[max][j])){
                max = k; // 得到 j 列的最大系数所在的行
            }
        }
        // 最大元系数调换到第 r 行
        swap(a[r], a[max]);
        if(fabs(a[r][j]) < eps) continue;
        // 主元系数变为 1
        for(int k = n + 1; k >= 1; k--){
            a[r][k] = a[r][k] / a[r][j];
        }
        // 其余行系数变成 0 
        for(int i = 1; i <= n; i++){
            if(i == r) continue;
            double tmp = a[i][j];
            for(int k = 1; k <= n + 1; k++){
                a[i][k] = a[i][k] - tmp * a[r][k];
            }
        }
        r ++;
    }
    r --;
    bool vis1 = false, vis2 = false;
    for(int i = 1; i <= n; i++){
        bool ok = false;
        for(int j = 1; j <= n; j++){
            if(fabs(a[i][j]) > eps) ok = true;
        }
        if(!ok && fabs(a[i][n + 1]) < eps) vis1 = true;
        if(!ok && fabs(a[i][n + 1]) > eps) vis2 = true;
    }
    if(vis2) return cout << "-1" << '\n', 0;
    if(vis1) return cout << "0" << '\n', 0;
    for(int i = 1; i <= n; i++){
        cout << "x" << i;
        printf("=%.2lf\n", a[i][n + 1]);
    }
    return 0;
}

1|3关于多项异或线性方程式求解

XOR Neighbors
观察每个点的性质, 发现对于任何一个点 $u$ 来说, 设其邻点为 $v$ , 则有 $⨁_{v_{i} \in Adj (u)} x_{v_{i}} = 0$ , 那么这样其实就符合一个异或方程组, 对于 $n$ 个点就有 $n$ 个方程组, 那么直接使用高斯消元求, 但是异或方程式与普通的方程式并不一样, 这里介绍一下异或方程主要特征:

在消除当前行的其余元时, 应逆序在这一行上异或上主元乘当前元, 这样可以使当前行的元变为 $0$
异或方程式无法像一般线性方程式那样可以化简成行最简形式, 经过高斯消元之后会变成上三角的样式
在求解时, 应从最后一行一步一步向上求解, 由于不是行最简, 所以要把所有的非零元代入方程求解

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
signed main()
{
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    int n, m; cin >> n >> m;
    vector<vector<int>> a(n + 10, vector<int> (n + 10));
    for(int i = 1; i <= m; i++){
        int u, v; cin >> u >> v;
        a[u][v] = a[v][u] = 1;
    }
    int r = 1;
    for(int j = 1; j <= n; j++){
        int mx = r;
        for(int i = r + 1; i <= n; i++){
            if(a[i][j] > a[mx][j]){
                mx = i;
            }
        }
        swap(a[r], a[mx]);
        for(int i = r + 1; i <= n; i++){
            if(i == r || a[i][j] == 0) continue;
            for(int k = n + 1; k >= j; k--){
                a[i][k] = a[i][k] ^ (a[i][j] * a[r][k]);
            }
        }
        r++;
    }
    int tot = 0;
    for(int i = n; i >= 1; i--){
        for(int j = i + 1; j <= n; j++){
            a[i][n + 1] = a[i][n + 1] ^ (a[i][j] * a[j][n + 1]);
        }
        if((a[i][i] == 0 && a[i][n + 1] != 0) || (a[i][i] != 0 && a[i][n + 1] == 0)){
            cout << "No" << '\n';
            return 0;
        }
        if(a[i][i] == 0 && a[i][n + 1] == 0) a[i][n + 1] = (1 << tot++);
    }
    cout << "Yes" << '\n';
    for(int i = 1; i <= n; i++) cout << a[i][n + 1] << ' ';
    return 0;
}

2|0例题

Painter's Problem

题目简述: 有一个由 n*n 个小方砖组成的正方形墙。一些砖块是白色的，而一些是黄色的。
鲍勃是一名画家，他想把所有的砖块都涂成黄色。但鲍勃的刷子有点问题，一旦他用这个刷子涂上砖块 (i, j)，
那么 (i, j)、(i-1, j)、(i+1, j)、(i, j-1) 和 (i, j+1) 的砖块颜色都会发生变化。
你的任务是找到鲍勃应该涂的最少砖块数量，以使所有的砖块都变成黄色。

本题乍一看很难想到使用高斯消元，下面通过建模转化成线性方程组问题。
共有 $k = n * n$ 个方砖，那么设 $x_{i}$ 为当前方砖刷的次数，明显的，每个方砖刷的次数最多为 $1$ 。

定义 $a_{i, j}$ 表示两块砖的关系，若刷了第 $j$ 块砖头对 $i$ 有影响那么就是 $1$ ，特别的 $a_{i i} = 1$ 。

对于方砖 $i$ 的颜色，设黄色为 $0$ ，白色为 $1$ ，那么最后颜色取决于初始颜色 $s_{i}$ 和所有其 $a_{i j} = 1$ 的异或操作情况

故有如下方程组。

{\begin{cases} s_{1} \oplus a_{1, 1} & & x_{1} \oplus a_{1, 2} & & x_{2} \oplus \dots \oplus a_{1, k} & & x_{k} = 0 \\ s_{2} \oplus a_{2, 1} & & x_{1} \oplus a_{2, 2} & & x_{2} \oplus \dots \oplus a_{2, k} & & x_{k} = 0 \\ \dots \\ s_{k} \oplus a_{k, 1} & & x_{1} \oplus a_{k, 2} & & x_{2} \oplus \dots \oplus a_{k, k} & & x_{k} = 0 \end{cases}

那么 $x_{i} = 1$ 就是被刷了一次，那么最后的答案就是 $\sum_{i = 1}^{k} x_{i}$

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string.h>
#include<algorithm>
#include<map>
#include<queue>
#include<vector>
#include<cmath>
#include<stdlib.h>
#include<time.h>
using namespace std;
const int N = 300, mod = 1e9 + 7;
int a[N][N], dx[] = {0, 0, 1, -1}, dy[] = {1, -1, 0, 0};
void solve(){
    int n; cin >> n;
    memset(a, 0, sizeof a);
    for(int i = 1; i <= n * n; i ++){
        char x; cin >> x;
        if(x == 'y') a[i][n * n + 1] = 0;
        else if(x == 'w') a[i][n * n + 1] = 1;
    }
    for(int i = 1; i <= n; i ++){
        for(int j = 1; j <= n; j++){
            a[n * (i - 1) + j][n * (i - 1) + j] = 1;
            for(int k = 0; k < 4; k ++){
                int xx = i + dx[k], yy = j + dy[k];
                if(xx < 1 || xx > n || yy < 1 || yy > n) continue;
                a[n * (i - 1) + j][n * (xx - 1) + yy] = 1;
            }
        }
    }
    int r = 1;
    for(int j = 1; j <= n * n; j ++){
        int max = r;
        for(int k = r + 1; k <= n * n; k ++){
            if(abs(a[k][j]) > abs(a[max][j])) max = k;
        }
        if(a[max][j] == 0) continue;
        for(int k = 1; k <= n * n + 1; k ++) swap(a[r][k], a[max][k]);
        for(int i = 1; i <= n * n; i ++){
            if(i == r || a[i][j] == 0) continue;
            for(int k = 1; k <= n * n + 1; k ++) a[i][k] ^= a[r][k];
        }
        r ++;
    }
    for(int i = r; i <= n * n; i++){
        if(a[i][n * n + 1]) return cout << "inf" << '\n', void();
    }
    int res = 0;
    vector<int> x(N);
    for(int i = 1; i <= n * n; i ++){
        x[i] = a[i][n * n + 1];
        for(int j = i + 1; j <= n * n; j++) x[i] ^= (a[i][j] && x[j]);
        res += x[i];
    }
    cout << res << '\n';
}
signed main()
{
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    int t; cin >> t; while(t--) solve();
    return 0;
}

开关问题
与上述的异或方程式同理, 不再赘述, 每个灯都能影响其他的灯, 所以可以列出对所有对该灯有影响的异或方程式: $⨁_{A_{i, j} * x_{i}} = e n d_{i}$ , 对于一个自由变量, 那么我们就有开或不开两种, 乘 $2$ 即可

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6 + 10, mod = 1e9 + 7;
void solve(){
    int n; cin >> n;
    vector<vector<int>> a(n + 10, vector<int> (n + 10));
    for(int i = 1; i <= n; i++) cin >> a[i][n + 1];
    for(int i = 1; i <= n; i++){
        int x; cin >> x;
        a[i][n + 1] ^= x;
        a[i][i] = 1;
    }
    int x, y;
    while(cin >> x >> y && (x || y)){
        a[y][x] = 1;
    }
    int r = 1;
    for(int j = 1; j <= n; j++){
        int mx = r;
        for(int i = r + 1; i <= n; i++){
            if(a[i][j] > a[mx][j]){
                mx = i;
            }
        }
        swap(a[r], a[mx]);
        for(int i = r + 1; i <= n; i++){
            if(i == r || a[i][j] == 0) continue;
            for(int k = n + 1; k >= j; k--){
                a[i][k] = a[i][k] ^ (a[i][j] * a[r][k]);
            }
        }   
        r++;
    }  
    int res = 1;
    for(int i = n; i >= 1; i--){
        for(int j = i + 1; j <= n; j++){
            a[i][n + 1] = a[i][n + 1] ^ (a[i][j] * a[j][n + 1]);
        }
        if(a[i][i] == 0 && a[i][n + 1] != 0){
            return cout << "Oh,it's impossible~!!" << '\n', void();
        }
        if(a[i][i] == 0 && a[i][n + 1] == 0) res *= 2;
    }
    cout << res << '\n';
}
signed main(){
    std::ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int t;cin>>t;while(t--)solve();
}

__EOF__

本文作者：Sakurajimamai
本文链接：https://www.cnblogs.com/o-Sakurajimamai-o/p/18324342.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！