[线性代数]高斯消元法

前言

仍然是自闭的一天。
学一学搞死高斯消元弄昏一下脑子。

概述

高斯消元主要用于求解线性方程组，也可以求矩阵的秩，矩阵的逆等。该算法的时间复杂度与方程组个数，未知数个数等有关。一般的时间复杂度为$O(n^3)$

线性方程组

含有多个一次未知数的方程组，形如
{ a 1 , 1 x 1 + a 1 , 2 x 2 + a 1 , 3 x 3 + . . . + a 1 , m x m = b 1 a 2 , 1 x 1 + a 2 , 2 x 2 + a 2 , 3 x 3 + . . . + a 2 , m x m = b 2 a 3 , 1 x 1 + a 3 , 2 x 2 + a 3 , 3 x 3 + . . . + a 3 , m x m = b 3 ⋯ a n , 1 x 1 + a n , 2 x 2 + a n , 3 x 3 + . . . + a n , m x m = b n \left\{ \begin{aligned} &a_{1,1}x_1+a_{1,2}x_2+a_{1,3}x_3+...+a_{1,m}x_m =b_1\\ &a_{2,1}x_1+a_{2,2}x_2+a_{2,3}x_3+...+a_{2,m}x_m =b_2\\ &a_{3,1}x_1+a_{3,2}x_2+a_{3,3}x_3+...+a_{3,m}x_m =b_3\\ & \quad \quad\quad\quad\quad\quad\quad\quad\quad \cdots \\ &a_{n,1}x_1+a_{n,2}x_2+a_{n,3}x_3+...+a_{n,m}x_m =b_n\\ \end{aligned} \right. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​​a1,1​x1​+a1,2​x2​+a1,3​x3​+...+a1,m​xm​=b1​a2,1​x1​+a2,2​x2​+a2,3​x3​+...+a2,m​xm​=b2​a3,1​x1​+a3,2​x2​+a3,3​x3​+...+a3,m​xm​=b3​⋯an,1​x1​+an,2​x2​+an,3​x3​+...+an,m​xm​=bn​​
记作矩阵形式：

$A=\left[\begin{array}{ccccc}{a}_{1,1}& a_{1,2}& a_{1,3}& \cdots & a_{1,m}\\ a_{2,1}& a_{2,2}& a_{2,3}& \cdots & a_{2,m}\\ a_{3,1}& a_{3,2}& a_{3,3}& \cdots & a_{3,m}\\ \cdots & \cdots & \cdots & \cdots & \cdots \\ a_{n,1}& a_{n,2}& a_{n,3}& \cdots & a_{n,m}\end{array}\right],B=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ \cdots \\ b_n\end{array}\right],Ax=B$

单位矩阵

在矩阵的乘法中，有一种矩阵起着特殊的作用，如同数的乘法中的1，这种矩阵被称为单位矩阵。它是个方阵，从左上角到右下角的对角线（称为主对角线）上的元素均为1，除此以外全都为0。

根据单位矩阵的特点，任何矩阵与单位矩阵相乘都等于本身。
单位矩阵

增广矩阵

又称扩增矩阵，就是在系数矩阵的右边添上一列，这一列是线性方程组的等号右边的值。

矩阵的初等变换

矩阵的初等变换分为矩阵的初等行变换和矩阵的初等列变换。矩阵的初等行变换和矩阵的初等列变换统称为初等变换。

等价

如果矩阵$B$可以由$A$通过一系列初等变换得到，则称矩阵$A$和$B$等价。

初等行变换

数域$P$上的初等行变换的定义：下面三种变换：

1. 以$P$中某个非$0$数乘矩阵的某一行
2. 把矩阵某一行的$C$倍加到另一行，$C$为P中任意一个数
3. 互换矩阵某两行的位置

若矩阵$A$通过初等行变换变成了矩阵$B$，写作$A\to B$

任意矩阵通过一系列初等行变换都能转化为阶梯矩阵；若该矩阵为方阵，则为上三角矩阵。

初等列变换

与初等行变换进行类比，得到定义：

1. 以$P$中某个非$0$数乘矩阵的某一列
2. 把矩阵某一行的$C$倍加到另一列，$C$为P中任意一个数
3. 互换矩阵某两列的位置

初等矩阵

定义：由单位矩阵$E$经过一系列初等变换后得到的矩阵$E^{\prime }$，称作初等矩阵。
引理：
对一个$n\ast m$的矩阵$A$作一次初等行变换，等价于在$A$左边乘上一个$n\ast n$的相应初等矩阵
对一个$n\ast m$的矩阵$A$作一次初等列变换，等价于在$A$右边乘上一个$m\ast m$的相应初等矩阵

高斯消元法

具体步骤

得到一个方程组。例如：
$$\{\begin{array}{rl}& 2x_1+1x_2+1x_3=1\\ & 4x_1+7x_2+8x_3=31\\ & 5x_1+2x_2+1x_3=4\end{array}$$
将它转换为增广矩阵：
$$\left[\begin{array}{cccc}2& 1& 1& 1\\ 4& 7& 8& 31\\ 5& 2& 1& 4\end{array}\right]$$

1. 首先，找到第一个第一列不为$0$的行(一般找绝对值最大的，如果都为$0$就跳过当前步)。
2. 然后把这个数带入下面的方程组，把下面的第一列换掉代入，那么下面方程组的第一列都化为$0$了。

例如这个方程，先找最大的第一列
$$\left[\begin{array}{cccc}5& 2& 1& 4\\ 4& 7& 8& 31\\ 2& 1& 1& 1\end{array}\right]$$
(具体的步骤就不说了)代入消元：
$$\left[\begin{array}{cccc}5& 2& 1& 4\\ 0& \frac{27}{5}& \frac{36}{5}& \frac{139}{5}\\ 0& \frac{1}{5}& \frac{3}{5}& -\frac{3}{5}\end{array}\right]$$
3. 除了最后一列，对于每列都这样找一行消掉。因为有解的方程组是个方阵，最终会转换为上三角矩阵。
消成上三角矩阵：
$$\left[\begin{array}{cccc}5& 2& 1& 4\\ 0& \frac{27}{5}& \frac{36}{5}& \frac{139}{5}\\ 0& 0& -\frac{9}{27}& -\frac{776}{135}\end{array}\right]$$
很显然，就可以从下往上，依次得到$x_3,x_2,x_1$了。

Code：

普通的(湘妹儿的)：

bool gauss(int n,db a[][M])
{
	bool flag=1;
	for (int i = 1;i <= n; ++ i)
	{
		int Max = i;
		for (int r = i + 1; r <= n; ++ r)
			if (cmp(a[Max][i], a[r][i]) == -1)
				Max = r;
		swap(a[i], a[Max]);
		if (cmp(a[i][i], 0) == 0)
		{
			flag = 0;
			continue;
		}
		for (int j = 1; j <= n; ++ j)
		{
			if (i == j || cmp(a[j][i], 0) == 0)
				continue;
			for (int k = i + 1; k <= n + 1; ++ k)
				a[j][k] -= a[i][k] * a[j][i] / a[i][i];
			a[j][i] = 0;
		}
	}
	return flag;
}

异或的：
```cpp
LL Gauss(LL N,LL M)
{
	LL tra, col;
    for (tra = 1, col = 1; tra <= N && col <= M; ++ tra, ++ col)
    {
        int Maxx = tra;
        for (Int i = tra; i <= N; ++ i)
            if (B[i][col])
            {
            	Maxx = i;
            	break;
			}
        if (! B[Maxx][col])
		{
            tra --;
            continue;
        }
        if (Maxx != tra)
            for (Int i = col; i <= M + 1; ++ i)
				Swap(B[Maxx][i], B[tra][i]);
        for (Int i = tra + 1; i <= N; ++ i)
		{
            if (! B[i][col])
				continue;
            for (Int j = col; j <= M + 1; ++ j)
                B[i][j] ^= B[tra][j]; 
        }
    }
    for (Int i = tra; i <= N; ++ i)
        if (B[i][col])
			return -1;
    return N - tra;
}

高斯-约旦消元法

特点

精度更高，代码更短，不用回带。

基本思想

模拟加减消元，将每个方程转化为$k_i{x}_i=b_i$的形式，然后带回求得每个$x_i$

具体思路

1. 选择一个没被解决过的未知数作为主元。再选择包含这个未知数的一个方程
2. 将这个方程主元的系数化为$1$
3. 将消元的式子带入其他方程加减消元，消去这个主元。(包括当前方程，这也是高斯-约旦消元法和普通高斯消元法的区别所在)
4. 重复步骤，将方程组化为对角矩阵(即系数在对角线上，而普通高斯消元则是转化为上三角矩阵)。然后求解即可。
   其余操作与普通高斯消元相似。

Code

题目

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define DB double
#define MAXN 105
#define LL long long
#define Int register int
using namespace std;
DB A[MAXN][MAXN];
DB FABS(DB x)
{
	if (x > 0)
		return x;
	return -x;
}
void Swap(DB &x,DB &y)
{
	DB temp = x;
	x = y;
	y = temp;
}
int main()
{
	int n;
	scanf("%d", &n);
	for (Int i = 1; i <= n; ++ i)
		for (Int j = 1; j <= n + 1; ++ j)
			scanf("%lf", &A[i][j]);
	for (Int i = 1; i <= n; ++ i) // 枚举列
	{
		int Maxx = i;
		for (Int j = i + 1; j <= n; ++ j) // 找绝对值最大的
			if (FABS(A[j][i]) > FABS(A[Maxx][i]))
				Maxx = j;
		for (Int j = 1; j <= n + 1; ++ j)
			Swap(A[Maxx][j], A[i][j]);
		if (! A[i][i]) // 没有这个未知数，无解（或无数解） 
			return !printf("No Solution\n");
		for (Int j = 1; j <= n; ++ j)
		{
			if (j != i)
			{
				DB temp = A[j][i] / A[i][i];
				for (Int k = i + 1; k <= n + 1; ++ k)
					A[j][k] -= A[i][k] * temp;
			}
		}
	}
	for (Int i = 1; i <= n; ++ i)
		printf("%.2lf\n", A[i][n + 1] / A[i][i]);
	return 0;
}