【OI学习笔记】线性代数-高斯消元法

宁可少些，但要好些，二分之一个证明等于 0
——德国数学家、物理学家 高斯

板块：线性代数-高斯消元法
难度：较难
前置知识点：矩阵的基本概念和三角阵、矩阵的初等变换
前置知识一览：

• 矩阵的概念和特殊矩阵
• 矩阵的初等变换（待补）

高斯消元法是线性代数中的重点知识，也是信息学奥赛中的考纲知识点，一般用于求解形如下所示的多元线性方程组。
$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{3n}{x}_n=b_2\\ ⋮\\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=b_n;\end{array}$$

高斯消元法的简单原理概述

概述： 高斯消元法的基本原理是把多元线性方程组的系数和等式右边记录为一个矩阵，并通过矩阵的变换使其成为上三角阵。最后再判定根据行数判断其有无解的情况，如果有解，还需要倒序回溯求解。

转换为上三角阵

• 【模板】高斯消元解线性方程组
  以该方程组为例：
  $$\{\begin{array}{l}{x}_1+2x^2-x_3=-6\\ 2x_1+x_2-3x_3=-9\\ -x_1-x_2+2x_3=7\end{array}$$
  通过手动计算可以得到该方程组的解为 $\{\begin{array}{l}{x}_1=1\\ x_2=-2\\ x_3=3\end{array}$.将其方程组左边的系数和等式右边转化为如下矩阵：
  $$\left(\begin{array}{cccc}1& 2& -1& -6\\ 2& 1& -3& -9\\ -1& -1& 2& 7\end{array}\right)$$
  高斯消元法是基于初等行变换的，也就是如下 3 条：
• 把某一行乘一个非零的数（操作①）
• 交换某两行（操作②）
• 把某行的若干倍加到另一行上（操作③）
  通过这三条规则的变换，最后应该要得到如下的上三角阵，但此时还没有求出解，因为各方程组还存在多个未知数：（具体变换过程将在下一项算法步骤中详细呈现）：
  $$\left(\begin{array}{cccc}1& 0.5& -1.5& -4.5\\ 0& 1& \frac{1}{3}& -1\\ 0& 0& \frac{2}{3}& 3\end{array}\right)$$

判定有无解

完美阶梯形：在执行完变换上三角阵的操作后，恰好为 $n$ 行且 $1\sim n$ 的各方程未知数个数为 $n,n-1,\cdots ,1$

• 如果恰好为完美阶梯形，则有唯一解。
• 如果不是完美阶梯形：
  • 如果出现了左边无未知数（为 0）但右边系数非 0 的矛盾情况，则说明该方程组无解。
  • 否则就是有 $0=0$ 的情况，则该方程组有 $n$ 个未知数， $n-1$ 个方程，有无穷组解。

比如对于上文给出的样例，就是有唯一解的情况。

计算答案

对于无解或有无穷组的情况，我们根据题目要求输出即可，但对于有唯一解的情况，我们必然是要求出唯一解的。

求解的过程就是顺次计算 $x_n,x_{n-1},x_{n-2},\cdots ,x_1$。因为我们已经将矩阵化为了上三角阵，所以显然，求解的过程是先计算 $x_n$ 所在的一元一次方程，然后将 $x_n$ 代入到上一行，就可以只计算 $x_{n-1}$ 了，再将 $x_n,x_{n-1}$ 代入 上一行，计算 $x_{n-2}\cdots$

最后，我们会将矩阵化为如下的形式，也就是除各行需要求解的未知数的系数外的其他系数均化为了 0：
$$\left(\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 0& -2\\ 0& 0& 1& 3\end{array}\right)$$
以上是高斯消元法的简单说明，接下来将结合算法步骤详细讲述全过程
。

算法步骤及详细过程

高斯消元法的算法步骤如下，在步骤 1 中，经过一轮操作后，会固定一个方程，最后变换成的上三角阵就是经过 $n$ 轮变换后的矩阵：

1. 从已经枚举每一列 $c$
   （1）找到第 $c$ 列绝对值最大的一行（为防止精度丢失）
   （2）将该行换至上一次固定的行的下一行（对应操作 ②）
   （3）将该行第 1 个数数变成 1（对应操作 ①）
   （4）将下面所有行的第 $c$ 列消为 0（对应操作 ③）
2. 判断有无解
3. 在有解的情况下自后向前求解（依次代入变换各方程为一元一次方程）

依然以该矩阵为例：
$$\left(\begin{array}{cccc}1& 2& -1& -6\\ 2& 1& -3& -9\\ -1& -1& 2& 7\end{array}\right)$$
从第 0 列开始枚举，定义 $r$ 表示上一次固定的行编号加 1，初始 $r=0$。然后开始执行上述算法步骤。

步骤 1 的（1）找到绝对值最大的一行，也就是第二行，将它换至第 0 行，也就是：
$$\left(\begin{array}{cccc}2& 1& -3& -9\\ 1& 2& -1& -6\\ -1& -1& 2& 7\end{array}\right)$$
然后，将该行的第 1 个数变成 1，也就是让每一个数都除以 2，变成：
$$\left(\begin{array}{cccc}1& 0.5& -1.5& -4.5\\ 1& 2& -1& -6\\ -1& -1& 2& 7\end{array}\right)$$
最后，将下面每一行的第 $c$ 列化为 0，那么就可以运用操作③，将第二行减去第一行，变成：
$$\left(\begin{array}{cccc}1& 0.5& -1.5& -4.5\\ 0& 1.5& 0.5& -1.5\\ 0& -0.5& 0.5& 2.5\end{array}\right)$$
然后再枚举第 $2\sim n-1$ 列，反复迭代，即可得到上文中给出的最终上三角阵：
$$\left(\begin{array}{cccc}1& 0.5& -1.5& -4.5\\ 0& 1& \frac{1}{3}& -1\\ 0& 0& \frac{2}{3}& 3\end{array}\right)$$
这一上三角阵恰好有 $n$ 行且是一个完美阶梯形，因此可以判断是具有唯一解的，按照讲解原理时的步骤计算，先将 $x_n$ 计算出来，将其系数化为 1，得：
$$\left(\begin{array}{cccc}1& 0.5& -1.5& -4.5\\ 0& 1& \frac{1}{3}& -1\\ 0& 0& 1& 3\end{array}\right)$$
得到 $x_3=3$，然后再向上回溯，第 2 行减去第 3 行乘 $\frac{1}{3}$，将第 2 行的 $x_3$ 的系数化为 0（消元），得到：
$$\left(\begin{array}{cccc}1& 0.5& -1.5& -4.5\\ 0& 1& 0& -2\\ 0& 0& 1& 3\end{array}\right)$$
不难发现，第 2 个方程也解出来了，得 $x_2=-2$。再考虑第一行，先消掉 0.5，将第一行加上第二行乘0.5；再消掉 -1.5，将第一行加上第三行乘 -1.5（因为 $a_{23}$ 为0，无法消元），得：
$$\left(\begin{array}{cccc}1& 0& 0& 1\\ 0& 1& 0& -2\\ 0& 0& 1& 3\end{array}\right)$$
于是，我们得到了该方程组的唯一解为 $\{\begin{array}{l}{x}_1=1\\ x_2=-2\\ x_3=3\end{array}$.

特别注意，在实际算法中，由于可能出现因为 C++ 的精度问题导致的变量存储不准确，比如 0 存为0.00000…01 这种情况，因此我们判定为 0 的情况，就是判定某个数是否小于一个很小的数，比如 $eps=1e-8$。

参考代码

有注释版：

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

const int N = 110;
const double eps = 1e-9;

int n;
double a[N][N];

int gauss()
{
    int c, r;
    for (c = 0, r = 0; c < n; c++)
    {
        int t = r; //上一次固定的行数+1
        for (int i = r; i < n; i++)
        {
        	//找到第c行绝对值最大的一列
            if (fabs(a[i][c]) > fabs(a[t][c]))
             t = i;
        }
        
       //如果该行的第c列为0就无意义了，后面会出现除以0的情况
        if (fabs(a[t][c]) < eps) continue;
        
        for (int i = c; i <= n; i++) swap(a[t][i], a[r][i]); //换行
        for (int i = n; i >= c; i--) a[r][i] /= a[r][c]; //将第一个数变成1
        //将下面的每一行的第c列都变成0
        for (int i = r + 1; i < n; i++)
        {
            if (fabs(a[i][c]) > eps)
            {
                for (int j = n; j >= c; j--)
                {
                    a[i][j] -= a[r][j] * a[i][c];
                }
            }
        }
        r++; //加1，下一轮操作时就可以直接用
    }
    
    //非完美矩形
    if (r < n)
    {
        for (int i = r; i < n; i++)
        {
        	//右边出现了非0的情况
            if (fabs(a[i][n]) > eps)
                return 2; //无解
        }
        return 1; //无穷组解
    }
	//有解就要求解，将除等式右边外的矩形转化为一个对角阵
    for (int i = n; i >= 0; i--)
    {
        for (int j = i + 1; j <= n; j++)
        {
            a[i][n] -= a[i][j] * a[j][n];
        }
    }
    return 0;
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n + 1; j++)
            cin >> a[i][j];

    int ans = gauss();

    if (ans == 0)
    {
        for (int i = 0; i < n; i++)
        {
        	//可能出现-0.0的情况
            if (fabs(a[i][n]) < eps) a[i][n] = 0;
            printf("%.2lf\n", a[i][n]);
        }
    }
    else if (ans == 1) printf("Infinite group solutions"); //有无穷组解
    else printf("No solution"); //无解
    return 0;
}

无注释版：

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

const int N = 110;
const double eps = 1e-9;

int n;
double a[N][N];

int gauss()
{
    int c, r;
    for (c = 0, r = 0; c < n; c++)
    {
        int t = r;
        for (int i = r; i < n; i++)
        {
            if (fabs(a[i][c]) > fabs(a[t][c]))
             t = i;
        }
        
        if (fabs(a[t][c]) < eps) continue;
        
        for (int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
        for (int i = n; i >= c; i--) a[r][i] /= a[r][c];
        for (int i = r + 1; i < n; i++)
        {
            if (fabs(a[i][c]) > eps)
            {
                for (int j = n; j >= c; j--)
                {
                    a[i][j] -= a[r][j] * a[i][c];
                }
            }
        }
        r++;
    }
    
    if (r < n)
    {
        for (int i = r; i < n; i++)
        {
            if (fabs(a[i][n]) > eps)
                return 2;
        }
        return 1;
    }
    for (int i = n; i >= 0; i--)
    {
        for (int j = i + 1; j <= n; j++)
        {
            a[i][n] -= a[i][j] * a[j][n];
        }
    }
    return 0;
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n + 1; j++)
            cin >> a[i][j];

    int ans = gauss();

    if (ans == 0)
    {
        for (int i = 0; i < n; i++)
        {
            if (fabs(a[i][n]) < eps) a[i][n] = 0;
            printf("%.2lf\n", a[i][n]);
        }
    }
    else if (ans == 1) printf("Infinite group solutions"); 
    else printf("No solution");
    return 0;
}

*如何求解异或线性方程组

异或线性方程组中，未知数的系数均为 0 或 1，做的运算不是四则运算而是异或运算，因此，整体框架不变，只是有三处变动：

• 变动处1：判断是否为 0 的地方可以采用非运算。
• 变动处2：步骤1（4），将第 $r\sim n-1$ 行的第 $c$ 列变为 0的操作可以直接异或 $a[r][j]$。
• 变动处3：步骤 3，有解情况下需要求解的操作，将 $a[i][n]$ 异或上 $a[i][j]$ 与 $a[j][n]$ 做与运算的结果。

#include <iostream>
#include <cstdio>

using namespace std;

const int N = 110;

int n;
int a[N][N];

int gauss()
{
    int c, r;
    for (c = 0, r = 0; c < n; c++)
    {
        int t = r;
        for (int i = r; i < n; i++)
        {
            if (a[i][c])
            {
                t = i;
                break;
            }
        }
        
        if (!a[t][c]) continue;
        for (int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
        for (int i = r + 1; i < n; i++)
        {
            if (a[i][c])
            {
                for (int j = n; j >= c; j--)
                    a[i][j] ^= a[r][j];
            }
        }
        r++;
    }
    
    if (r < n)
    {
        for (int i = r; i < n; i++)
        {
            if (a[i][n])
                return 2;
        }
        return 1;
    }
    
    for (int i = n; i >= 0; i--)
    {
        for (int j = i + 1; j <= n; j++)
        {
            a[i][n] ^= a[i][j] & a[j][n];
        }
    }
    return 0;
}

int main()
{
    cin >> n;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n + 1; j++)
            cin >> a[i][j];
    
    int ans = gauss();
    
    if (ans == 0)
    {
        for (int i = 0; i < n; i++) printf("%d\n", a[i][n]);
    }
    else if (ans == 1) printf("Multiple sets of solutions");
    else printf("No solution");
    return 0;
}