浅谈高斯消元法

2023.6.16：发布

2023.8.29：修缮，加上自己觉得通俗易懂的理解，更新矩阵求逆。

高斯消元

高斯消元可以用于线性方程组求解或者行列式计算，求矩阵的逆等等，也算是比较基础的一个思想。

消元法

定义

消元法是将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其带入到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的。消元法主要用于二元一次方程组的求解。

举例：利用消元法求解二元一次线性方程组：

\[\left\{\begin{matrix} 3x+2y=110\\ x-y=10 \end{matrix}\right. \]

我们很容易由二式得到：

\[x=y+10 \]

我们把这个式子代入一式，然后得到：

\[3(y+10)+2y=110 \]

展开得到：

\[5y+30=110 \]

然后一式的未知数由两个变为了一个，我们可以解出：

\[y=16 \]

然后代入一式或者二式就可以解出来 \(x\) 的值了。

消元法理论核心

• 两方程互换，解不变

• 一方程乘以非零数 \(k\)，解不变

• 一方程乘以数 \(k\) 加上另一方程，解不变（此时的 \(k\) 可以是 \(0\)）

高斯消元法思想概念

德国数学家高斯对消元法进行了思考分析，得出了如下结论：

• 在消元法中，参与计算和发生改变的是方程中各变量的系数；

• 各变量并未参与计算，且没有发生改变；

• 可以利用系数的位置表示变量，从而省略变量；

• 在计算中将变量简化省略，方程的解不变。

高斯在这些结论的基础上，提出了高斯消元法，首先将方程的增广矩阵利用行初等变换化为行最简形，然后以线性无关为准则对自由未知量赋值，最后列出表达方程组通解。

高斯消元五步法

1. 增广矩阵行初等变换为行最简形

2. 还原线性方程组

3. 求解第一个变量

4. 补充自由未知量

5. 列表示方程组的通解

过程

我看到 OI Wiki 上写的太过复杂，所以我尽量用通俗易懂的语言来叙述。

例如我们需要解下面的线性方程组：

\[\left\{\begin{matrix} 2x_{1}+5x_{3}+6x_{4}=9\\ x_{3}+x_{4}=-4\\ 2x_{3}+2x_{4}=-8 \end{matrix}\right. \]

增广矩阵初等变换为行最简形

最简形矩阵

定义形如下面矩阵 \(A\) 的矩阵被称为最简形矩阵。

\[A=\begin{bmatrix} 1&0&5&6&9\\ 0&0&1&1&-4\\ 0&0&0&1&-8 \end{bmatrix} \]

其判断的依据可以简记为是否能画出一条阶梯状的线将矩阵分为上下两部分，满足在其下面的元素都是 \(0\)，特别的，画出的线的上方的元素只能是 \(1\) 或 \(0\)。

增广矩阵是由方程组系数矩阵 \(A\) 与常数列 \(b\) 的并生成的新矩阵，即 \((A|b)\)，增广矩阵行初等变换化为最简形，省略了变量而用变量的系数位置来表示变量，增广矩阵中用竖线隔开了系数矩阵和常数列，相当于等于号。

上面的线性方程组可以变换为下面的形式。

\[\begin{bmatrix} \begin{array}{cccc|c} 2&0&5&6&9\\ 0&0&1&1&-4\\ 0&0&2&2&-8 \end{array} \end{bmatrix} \]

我们用第三行减去两倍的第二行，得到：

\[\begin{bmatrix} \begin{array}{cccc|c} 2&0&5&6&9\\ 0&0&1&1&-4\\ 0&0&0&0&0 \end{array} \end{bmatrix} \]

然后第一行除以 \(2\) 得到：

\[\begin{bmatrix} \begin{array}{cccc|c} 1&0&2.5&3&4.5\\ 0&0&1&1&-4\\ 0&0&0&0&0 \end{array} \end{bmatrix} \]

然后我们用第一行减去 \(2.5\) 倍的第二行

\[\begin{bmatrix} \begin{array}{cccc|c} 1&0&0&0.5&14.5\\ 0&0&1&1&-4\\ 0&0&0&0&0 \end{array} \end{bmatrix} \]

至此化为了最简形矩阵

还原线性方程组

\[\left\{\begin{matrix} x_{1}+0.5x_{4}=14.5\\ x_{3}+x_{4}=-4 \end{matrix}\right. \]

也就是把最简形矩阵重新书写成线性方程组的形式。

求解第一个变量

\[\left\{\begin{matrix} x_{1}=-0.5x_{4}+14.5\\ x_{3}=-x_{4}-4 \end{matrix}\right. \]

也就是将每一行第一个系数不为 \(0\) 的变量用其他变量表达出来，也就是使左边只有需要表达的变量且系数为 \(1\) ，然后其他的项移到右边。

补充自由未知量

\[\left\{\begin{matrix} x_{1}=-0.5x_{4}+14.5\\ x_{2}=x_{2}\\ x_{3}=-x_{4}-4\\ x_{4}=x_{4} \end{matrix}\right. \]

第三步中得到的方程式只有 \(x_{1}\) 和 \(x_{3}\) 的，说明其他的变量都是不受其他变量约束的，也就是可以取任意值。所以只能是自己等于自己。

列表示方程组的通解

\[\begin{bmatrix} x_{1}\\x_{2}\\x_{3}\\x_{4} \end{bmatrix}= \begin{bmatrix} 0\\1\\0\\0 \end{bmatrix}x_{2}+ \begin{bmatrix} -0.5\\0\\-1\\1 \end{bmatrix}x_{4}+ \begin{bmatrix} 14.5\\0\\-4\\0 \end{bmatrix}= \begin{bmatrix} 0\\1\\0\\0 \end{bmatrix}c_{1}+ \begin{bmatrix} -0.5\\0\\-1\\1 \end{bmatrix}c_{2}+ \begin{bmatrix} 14.5\\0\\-4\\0 \end{bmatrix} \]

其中的 \(c_{1},c_{2}\) 为任意常数。

由于 \(x_{2}\) 和 \(x_{4}\) 是自由未知量，所以可以随便取值，所以在解的右边令二者分别为任意常数即完成对方程组的求解。

其实可以看做化为一个上三角矩阵然后从下往上依次回代的过程。

我们就是使每一行的对角线上的元素都是 \(1\) 这样刚好对应了每一个变量为 \(1\) 到了最后一行的时候，会获得一个变量的值，然后依次往上回代即可。最后即可求出所有变量的值。

在处理矩阵的时候，我们一般用到以下几种操作：

• 交换两行

• 系数化为 \(1\)

• 某行加上或减去另一行的 \(k\) 倍

首先我们确定第 \(i\) 行选定 \(x_{i}\) 为主元，然后往下找第 \(i\) 列不是零的数，然后交换这两行，将交换后的第 \(i\) 行主元系数化为 \(1\)

然后往下遍历，下面行只要第 \(i\) 列还有值，就加上减去第 \(i\) 行的 \(k\) 倍使其系数化为 \(0\)，如此反复即可。

整个高斯消元的过程：首先我们钦定每一个变量系数为 \(1\) 的行在第 \(i\) 行，我们找到一行 \(x_{i}\) 系数不为 \(0\) 的一行，然后再枚举把下面每一行的所有 \(x_{i}\) 的系数都化为 \(0\)，然后到了最后得到了一个上三角矩阵，这样我们从最后一行，依次把值往上面的方程组中回代，逐步求得所有方程的解。

P3389【模板】高斯消元法 - 洛谷

参考代码：

/*
 * @Author: Aisaka_Taiga
 * @Date: 2023-08-29 11:02:40
 * @LastEditTime: 2023-08-29 16:48:06
 * @LastEditors: Aisaka_Taiga
 * @FilePath: \Desktop\P3389.cpp
 * 心比天高，命比纸薄。
 */
#include<bits/stdc++.h>

#define DB double
#define N 1100

using namespace std;

int n;
DB a[N][N], x[N];

inline int work()
{
    for(int i = 1; i <= n; i ++)
    {
        int r = i;//当前变量
        for(int k = i; k <= n; k ++)//找到此列下面第一个不为0的值所在的行
            if(fabs(a[k][i]) > 0){r = k; break;}
        if(r != i) swap(a[r], a[i]);//交换行
        if(fabs(a[i][i]) == 0)return 1;//如果当前变量的系数是0则此变量可以取任意值，无数个解
        for(int j = n + 1; j >= i; j --)//将第i个变量所在的方程变成i变量系数为1的方程
            a[i][j] /= a[i][i];
        for(int k = i + 1; k <= n; k ++)//遍历后面的每一行 
            for(int j = n + 1; j >= i; j --)//遍历后面的每一列 
                a[k][j] -= a[k][i] * a[i][j];//把后面所有的第i个变量系数都变成0 
    }
    for(int i = n - 1; i; i --)//遍历每一行
        for(int j = i + 1; j <= n; j ++)//遍历每一列 
            a[i][n + 1] -= a[i][j] * a[j][n + 1];//a[i][j]是系数，a[j][n+1]是当前第j列变量的值，因为j-n的所有变量都算出来了所以直接用
    return 0;//有解
}

signed main()
{
    cin >> n;
    for(int i = 1; i <= n; i ++)
        for(int j = 1; j <= n + 1; j ++)
            cin >> a[i][j];
    int flag = work();
    if(flag != 0) return cout << "No Solution" << endl, 0;
    for(int i = 1; i <= n; i ++)
        printf("%.2lf\n", a[i][n + 1]);
    return 0;
}

高斯约旦消元法

他和普通的高斯消元的区别就是这个最后得到的矩阵是只有对角线是有值的，其余元素为 \(0\)。

最后得到的矩阵就像下面这样：

\[\begin{bmatrix} \begin{array}{ccc|c} a_{1}&\dots&0&c_{1}\\ \vdots&\ddots&\vdots&\vdots\\ 0&\dots&a_{n}&c_{n} \end{array} \end{bmatrix} \]

其中 \(c_{1}\) 到 \(c_{n}\) 为常数。

这样做的好处就是让每一个变量都直接乘上自己的系数得到一个值，这样就便于计算了。

最后直接循环一下计算即可。

当然弊端就是无法判断无解是没有解还是无数个解。

高斯约旦消元法和朴素的高斯消元法最大的区别就是最后得到的矩阵的不同，朴素的高斯消元法得到的就是一个上三角矩阵，我们想要得到具体的值还需要倒着一个一个往回带，但是高斯约旦最后得到的是一个对角矩阵，且每一个方程对应一个未知量，系数都为 \(1\)，不需要回代。

具体的过程就是，还是钦定第 \(i\) 行是 \(x_{i}\) 系数为 \(1\) 的方程，我们在选到这个方程之后，就把除第 \(i\) 行的 \(x_{i}\) 的系数都化为 \(0\)，这还是跟上面差不多的操作，由于前面的其他行的只有对应的变量系数为 \(1\)，而在其他行这个变量系数都是 \(0\)，所以不用担心系数会不是 \(1\)。搞完最后一行就得到了一个系数为单位矩阵的一堆方程组，而这个新方程组就是原来方程组的解。

参考代码

/*
 * @Author: Aisaka_Taiga
 * @Date: 2023-08-29 16:22:52
 * @LastEditTime: 2023-08-29 17:10:20
 * @LastEditors: Aisaka_Taiga
 * @FilePath: \Desktop\P3389高斯约旦消元法.cpp
 * 心比天高，命比纸薄。
 */
#include<bits/stdc++.h>

#define DB double
#define N 1100

using namespace std;

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

signed main()
{
	cin >> n;
	for(int i = 1; i <= n; i ++)
	  	for(int j = 1; j <= n + 1; j ++)
	    	cin >> a[i][j];
	for(int i = 1; i <= n; i ++)//? 枚举每一个变量
	{
		int pl = i;//? 当前行的变量编号
		for(int j = i; j <= n; j ++)//? 枚举每一行
		  	if(a[j][i] > a[pl][i]) pl = j;//?找到系数最大的一行
		if(a[pl][i] == 0) return cout << "No Solution" << endl, 0;//?如果是零就是无解
		swap(a[i], a[pl]);//? 交换两列
		DB k = a[i][i];//? 取出当前变量所在行的系数
		for(int j = 1; j <= n + 1; j ++)//? 枚举当前行的系数
		  	a[i][j] = a[i][j] / k;//? 系数化为1
		for(int j = 1; j <= n; j ++)//? 枚举每一行
		{
			if(i == j) continue;//? 如果是当前变量所在行就跳过，因为我们处理过了
			DB kk = a[j][i];//? 取出当前行的第i个变量的系数
			for(int o = i; o <= n + 1; o ++)//? 枚举当前行的每一个元素的系数
			    a[j][o] = a[j][o] - kk * a[i][o];//? 将除第i行以外的所有的i的行都化为1
		}
	}
	for(int i = 1; i <= n; i ++)
	 	printf("%.2lf\n", a[i][n + 1]);
	return 0;
}

高斯消元法对矩阵求逆

对于一个方阵 \(A\)，若存在一个矩阵 \(A'\) 使得 \(A\times A' =I\)，则这个方阵 \(A'\) 就是 \(A\) 的逆矩阵。

求解方法如下：

• 构造一个 \(n\times 2n\) 的矩阵，左半边是 \(A\)，右半边是 \(I_{A}\)。

• 然后对 \(n\times n\) 的左半部分进行高斯消元，最后的右半边矩阵就成了 \(A'\)。

要注意的是，如果左半部分最后得到的不是单位矩阵就说明不存在 \(A'\)。

下面来看一道模板题：

P4783 【模板】矩阵求逆

/*
 * @Author: Aisaka_Taiga
 * @Date: 2023-08-29 17:34:42
 * @LastEditTime: 2023-08-29 19:17:43
 * @LastEditors: Aisaka_Taiga
 * @FilePath: \Desktop\P4783.cpp
 * 心比天高，命比纸薄。
 */
#include <bits/stdc++.h>

#define int long long
#define P 1000000007
#define DB double
#define N 1010

using namespace std;

int n, a[N][N << 1];

inline int ksm(int a, int b)
{
    int res = 1;
    while(b)
    {
        if(b & 1) res = (res * a) % P;
        a = (a * a) % P;
        b >>= 1;
    }
    return res % P;
}

inline int Aisaka_Taiga()
{
    for(int i = 1; i <= n; i ++)
    {
        int p = i;
        for(int j = i + 1; j <= n; j ++)
            if(a[j][i] > a[p][i]) p = j;
        if(a[p][i] == 0) return 1;
        if(p != i) swap(a[p], a[i]);
        int k = ksm(a[i][i], P - 2);
        for(int j = 1; j <= 2 * n; j ++)
            a[i][j] = (a[i][j] * k) % P;
        for(int j = 1; j <= n; j ++)
        {
            if(i == j) continue;
            int k = a[j][i];
            for(int o = i; o <= 2 * n; o ++)
                a[j][o] = (a[j][o] - k * a[i][o] % P + P) % P;
        }
    }
    return 0;
}

signed main()
{
    cin >> n;
    for(int i = 1; i <= n; i ++)
        for(int j = 1; j <= n; j ++)
            cin >> a[i][j];
    for(int i = 1; i <= n; i ++) a[i][n + i] = 1;
    int ff = Aisaka_Taiga();
    if(ff == 1) return cout << "No Solution" << endl, 0;
    for(int i = 1; i <= n; i ++)
    {
        for(int j = n + 1; j <= 2 * n; j ++)
            cout << a[i][j] << " ";
        cout << endl;
    }
    return 0;
}