数代数

高斯消元与线性方程组

字数统计: 2.8k字   |   阅读时长: 13分
,

摘要: 本文介绍高斯消元的算法原理和代码模板。

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

本文介绍高斯消元的算法原理和代码模板。首先简要介绍一下线性方程组与高斯消元，以及高斯消元的特殊情况：无解和无穷解。最后通过两道题来看一下高斯消元的代码模板。

高斯消元与线性方程组

高斯消元是一种求解线性方程组的方法。假设线性方程组由 M 个 N 元一次方程共同组成，所有系数可以写成 M 行 N 列的系数矩阵。再加上常数列，可以写成 M 行 N + 1 列的增广矩阵。

解这个线性方程组的步骤可以概括成对增广矩阵的三类操作 (初等行变换)：

  1. 用一个非零数乘以某一行。
  2. 把其中一行的若干倍加到另一行
  3. 交换两行的位置

增广矩阵经过初等行变换最后得到阶梯型矩阵。阶梯型矩阵还可以再化简为简化阶梯型矩阵。

把增广矩阵变为简化阶梯型矩阵的线性方程组求解算法就是高斯消元算法。

核心思想是对每个未知量 $x_{i}$，找到一个 $x_{i}$ 系数非零，而 $x_{1},x_{2},…,x_{i-1}$ 的系数均为零的方程，然后用初等行变换把其它方程的 $x_{i}$ 的系数全都消成零。

高斯消元的特殊情况

  1. 高斯消元过程中可能会出现 0 = d 这种方程，d 是非零常数。方程组无解。
  2. 找不到一个 $x_{i}$ 系数非零，而 $x_{1},x_{2},…,x_{i-1}$ 的系数均为零的方程。方程组有无穷多组解，解由主元和自由元组成。下面是一个例子：
$$ \begin{gathered} \begin{bmatrix} 1 & 2 & -1 & 3 \\ 2 & 4 & -8 & 0 \\ -1 & -2 & 6 & 2 \\ \end{bmatrix} \Rightarrow \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 0 & -6 & -6 \\ 0 & 0 & 5 & 5 \\ \end{bmatrix} \Rightarrow \begin{bmatrix} 1 & 2 & -1 & 3 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix} \Rightarrow \begin{bmatrix} 1 & 2 & 0 & 4 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix} \end{gathered} $$

总结一下，在高斯消元完算法成后，判断方程组解的情况的方式如下：

  • 若存在系数全都为零，常数不为零的行，则方程组无解。
  • 若系数不全为零的行恰好有 N 个，说明主元有 N 个，方程组有唯一解。
  • 若系数不全为零的行有 k < N 个，说明主元有 k 个，自由元有 N - k 个，方程组有无穷多组解。

题目

207. 球形空间产生器

有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。

现在，你被困在了这个n维球体中，你只知道球面上n+1个点的坐标，你需要以最快的速度确定这个n维球体的球心坐标，以便于摧毁这个球形空间产生器。

注意： 数据保证有唯一解。

输入格式

1
2
3
4
5
第一行是一个整数n。

接下来的n+1行，每行有n个实数，表示球面上一点的n维坐标。

每一个实数精确到小数点后6位，且其绝对值都不超过20000。

算法：高斯消元

共 n + 1 个点，第 i 个点的 n 维坐标为 a[i][1],a[i][2],...,a[i][n], i = 1,2,…,n+1

求一个点 x，其 n 维度坐标为 x[1],x[2],...,x[n]，需要满足对每个 i = 1,2,…,n+1，都有 sum((a[i][j] - x[j]) ^ 2) = C C 为常数。该方程有 n + 1 个 n 元二次方程组成，不是线性方程组。

可以通过相邻两个方程做差，变成 n 个一元一次方程，同时消去常数 C。

写成标准的线性方程组

写出增广矩阵

$$ \begin{bmatrix} 2(a[1][1] - a[2][1]) & 2(a[1][2] - a[2][2]) & \cdots & 2(a[1][n] - a[2][n]) & \sum\limits_{j=1}\limits^{n}(a[1][j]^{2} - a[2][j]^{2})\\ 2(a[2][1] - a[3][1]) & 2(a[2][2] - a[3][2]) & \cdots & 2(a[2][n] - a[3][n]) & \sum\limits_{j=1}\limits^{n}(a[2][j]^{2} - a[3][j]^{2})\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 2(a[n][1] - a[n+1][1]) & 2(a[n][2] - a[n+1][2]) & \cdots & 2(a[n][n] - a[n+1][n]) & \sum\limits_{j=1}\limits^{n}(a[n][j]^{2} - a[n+1][j]^{2})\\ \end{bmatrix} $$

代码 (模板：保证有唯一解时的高斯消元)

$O(N^{3})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#include <vector>
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

const double EPS = 1e-8;

int main()
{
    int n;
    cin >> n;
    vector<vector<double>> a(n + 2, vector<double>(n + 1, 0.0));
    for(int i = 1; i <= n + 1; ++i)
        for(int j = 1; j <= n; ++j)
            cin >> a[i][j];

    vector<vector<double>> c(n + 1, vector<double>(n + 1, 0.0)); // 系数矩阵
    vector<double> b(n + 1, 0.0); // 常数列
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= n; ++j)
        {
            c[i][j] = 2 * (a[i][j] - a[i+1][j]);
            b[i] += a[i][j] * a[i][j] - a[i + 1][j] * a[i + 1][j];
        }

    // 高斯消元, 保证有唯一解
    for(int i = 1; i <= n; ++i)
    {
        // 找到 x[i] 的系数不为零的一个方程
        for(int j = i; j <= n ;++j)
        {
            if(fabs(c[j][i]) > EPS)
            {
                for(int k = 1; k <= n; ++k)
                    swap(c[i][k], c[j][k]);
                swap(b[i], b[j]);
                break;
            }
        }
        // 消去其它方程的 x[i] 的系数
        for(int j = 1; j <= n; ++j)
        {
            if(i == j) continue;
            double rate = c[j][i] / c[i][i];
            for(int k = i; k <= n; ++k)
                c[j][k] -= c[i][k] * rate;
            b[j] -= b[i] * rate;
        }
    }

    cout << std::fixed << std::setprecision(3);
    for(int i = 1; i <= n; ++i)
        cout << b[i] / c[i][i] << " ";
    cout << endl;
}

208. 开关问题

有N个相同的开关，每个开关都与某些开关有着联系，每当你打开或者关闭某个开关的时候，其他的与此开关相关联的开关也会相应地发生变化，即这些相联系的开关的状态如果原来为开就变为关，如果为关就变为开。

你的目标是经过若干次开关操作后使得最后N个开关达到一个特定的状态。

对于任意一个开关，最多只能进行一次开关操作。

你的任务是，计算有多少种可以达到指定状态的方法。（不计开关操作的顺序）

算法：高斯消元

x[i] 表示第 i 个开关的操作情况，x[i] = 1 表示按了这个开关，x[i] = 0 表示没有按。

a[i][j] 表示第 i 个开关与第 j 个开关的关联情况，a[i][j] = 1 表示按下 j 会影响 i、a[i][j] = 0 表示无影响，a[i][i] = 1

一个开关 i，其最终状态 dst[i]，取决于其最初状态 src[i]，以及所有与 i 有联系的开关的操作情况的异或和。

异或方程组

异或运算是不进位的加法，因此仍然可以写增广矩阵，其每个值为 0 或 1。高斯消元过程中将加法，减法改为异或，不需要执行乘法。

代码 (模板：高斯消元的特殊情况讨论)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include <iostream>
#include <vector>

using namespace std;

int main()
{
    int K;
    cin >> K;
    while(K--)
    {
        int N;
        cin >> N;

        vector<int> src(N + 1);
        for(int i = 1; i <= N; ++i)
            cin >> src[i];
        vector<int> dst(N + 1);
        for(int i = 1; i <= N; ++i)
            cin >> dst[i];

        // 系数矩阵
        vector<vector<int>> a(N + 1, vector<int>(N + 1, 0));
        int I, J;
        do{
            cin >> I;
            cin >> J;
            a[J][I] = 1;
        }while(I != 0 && J != 0);
        for(int i = 1; i <= N; ++i)
            a[i][i] = 1;

        // 常数列
        vector<int> b(N + 1);
        for(int i = 1; i <= N; ++i)
            b[i] = src[i] ^ dst[i];

        // 高斯消元
        int free_v = 0;
        for(int i = 1; i <= N; ++i)
        {
            // 第 i 个方程，第 i 个主元
            // 找 x[i] 系数不为零的一个方程
            bool get = false;
            for(int j = i; j <= N; ++j)
            {
                if(a[j][i] > 0) // 找到，将该方程交换到第 i 个方程的位置
                {
                    for(int k = 1; k <= N; ++k)
                        swap(a[i][k], a[j][k]);
                    swap(b[i], b[j]);
                    get = true;
                    break;
                }
            }

            // 判断系数是否全都为 0
            bool flag = true;
            for(int j = 1; j <= N; ++j)
                if(a[i][j] != 0)
                {
                    flag = false;
                    break;
                }
            if(flag && b[i] != 0) // 判断是否出现 0 = 1 的方程
            {
                // 出现 0 = 1 方程，无解
                free_v = -1;
                break;
            }

            // 判断 x[i] 是否为自由元
            if(!get)
            {
                ++free_v;
                continue;
            }

            // 消去其它方程 x[i] 的系数
            for(int j = 1; j <= N; ++j)
            {
                if(i == j) continue;
                if(a[j][i] == 1)
                {
                    for(int k = i; k <= N; ++k)
                        a[j][k] ^= a[i][k];
                    b[j] ^= b[i];
                }
            }
        }
        if(free_v == -1)
            cout << "Oh,it's impossible~!!" << endl;
        else if(free_v == 0)
            cout << 1 << endl;
        else
            cout << (1 << free_v) << endl;
    }
}

代码2 (模板：高斯消元将行向量压缩为一个整数)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <iostream>
#include <vector>

using namespace std;

int main()
{
    int K;
    cin >> K;
    while(K--)
    {
        int N;
        cin >> N;

        // 状态压缩后的增广矩阵
        vector<int> a(N + 1);
        for(int i = 1; i <= N; ++i)
            cin >> a[i];
        for(int i = 1; i <= N; ++i)
        {
            int x;
            cin >> x;
            a[i] ^= x;
        }

        int I, J;
        do{
            cin >> I;
            cin >> J;
            a[J] |= 1 << I; // a[J][I] = 1
        }while(I != 0 && J != 0);
        for(int i = 1; i <= N; ++i)
            a[i] |= 1 << i; // a[i][i] = 1

        // 高斯消元
        int ans = 1;
        for(int i = 1; i <= N; ++i)
        {
            // 第 i 个元 x[i]
            for(int j = i + 1; j <= N; ++j)
            {
                if(a[j] > a[i])
                    swap(a[i], a[j]);
            }

            // 判断系数是否全都为 0 : a[i] <= 1
            if(a[i] == 1) // 判断是否出现 0 = 1 的方程
            {
                // 出现 0 = 1 方程，无解
                ans = 0;
                break;
            }
            if(a[i] == 0)
            {
                // 出现 0 = 0 的方程
                ans = 1 << (N - i + 1);
                break;
            }

            // a[i] 最高位的 1 为主元，消去其它方程该位的系数
            for(int k = N; k >= 1; --k)
            {
                if(a[i] >> k & 1)
                {
                    for(int j = 1; j <= N; ++j)
                    {
                        if(i == j) continue;
                        if(a[j] >> k & 1)
                            a[j] ^= a[i];
                    }
                    break;
                }
            }
        }
        if(ans == 0)
            cout << "Oh,it's impossible~!!" << endl;
        else
            cout << ans << endl;
    }
}
Share