高斯消元 学习笔记

简介

高斯消元通常用于求解线性方程组。时间复杂度 \(\mathcal{O}(n^3)\)。

概念与核心

首先是一些很显然的性质：

1. 两方程互换，解不变。
2. 两方程相加，解不变。
3. 一方程乘以一非零实数 \(x\)，解不变。

在初中，我们是如何进行消元的？大多都是进行加减消元。如等式 \(x+2y=1\) 和 \(-x+y=2\)，我们会把它们相加，进而得到 \(3y=3\) 后进行求解。多元线性方程亦如此。

考虑我们要求解一下方程组。

\[\begin{cases}2x_1+9x_2-5x_3=10\\4x_1+20x_2+x_3=24\\x_1-\dfrac{1}{2}x_2+x_3=8\end{cases} \]

为了方便，我们采用一个矩阵的方式存下这个方程组的系数以及常数。

\[\begin{bmatrix}2& 9& -5 &|&10\\4&20&1 &|&24\\1&-\dfrac{1}{2}&1&|&8\end{bmatrix} \]

我们把这个矩阵叫做增广矩阵。

如果我们通过某些操作（初等行列变换），把这个增广矩阵化为一个阶梯状矩阵：

\[\begin{bmatrix}1& 5& \dfrac{1}{4} &|&6\\0&1&\dfrac{11}{2} &|&2\\0&0&1&|&\dfrac{13}{31}\end{bmatrix} \]

最后一行的解 \(x_3\) 我们就知道了，然后回带到倒数第二个方程里，那么 \(x_2\) 就知道了。最后回带到第一个方程里，那么 \(x_1\) 就解出来了。

也就是只要我们求出这么一个阶梯转矩阵，那么解就可以通过回带得出。

过程

我们一列一列的消元。假设我们当前消到了第 \(c\) 列，还有 \(r\sim n\) 行还需变为阶梯状矩阵（确实按理来说 \(r=c\)，但是会有特殊情况）。

1. 在该列选择一个绝对值最大的数，作为“主方程”。(为啥是绝对值最大？因为怕选到了 \(0\))
2. 把该方程与第 \(r\) 行的方程互换位置。（方便计算）
3. 把 \(r\) 行的第一个系数化为 \(1\)。（让矩阵的对角都为 \(1\)）
4. 把 \(r+1\sim n\) 行的第 \(c\) 列的系数都消为 \(0\)。（矩阵的左下角都为 \(0\)）。

重复这个过程。

具体地，在代码中，使用 \(g_{i,j}\) 记录增广矩阵的系数。

选择主方程：

int t=r;
	for(int i=r+1;i<=n;i++)
		if(fabs(g[i][c])>fabs(g[t][c]))
			t=i;
	if(fabs(g[t][c])<eps) continue;
	//如果选到了 0，那么可能是无数解，也有可能无解。
	//选到 0 后，直接跳过该列。

互换位置：

if(t!=r) swap(g[t],g[r]);

系数化 \(1\)：

for(int i=n+1;i>=c;i--) g[r][i]/=g[r][c]; 

后面的行该列系数消为 \(0\)：

for(int i=r+1;i<=n;i++)
	for(int j=n+1;j>=c;j--)
		g[i][j]-=g[r][j]*g[i][c];

为啥是这样？

代码中的 \(i\) 是用来枚举后面的行，\(j\) 是用来枚举每行的后面列。

第 \(r\) 行该列的值是 g[r][j]，减去的倍数是 -g[i][c]（因为 g[r][c]=1）。就得到代码中的式子。倒序是因为避免 g[i][c] 被过早更新，我们还需要使用它算倍数！

最后回带也是没有任何技术难度。

判定解

如果最后解出来是一个完美的阶梯状矩形，那么有唯一解。这是好想的。

关键是：如果它不是一个完美的阶梯状矩阵呢？比如：

\[\begin{bmatrix}1& 5& 4 &4 &|&6\\0&1&11 &3&|&2\\0&0&1&3&|&3\\ 0&0&0&0&|&1\\\end{bmatrix} \]

容易发现 \(0x=1\) 是无解的，也就是，算到最后 \(r\sim n\) 全是 \(0\)，如果右侧 \(n+1\) 列的常数有非 \(0\) 的，那么就是无解。

如果右侧全为 \(0\) 呢？\(0x=0\) 显然有无数解。

题目

P3389 【模板】高斯消元法

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n;
const int N=103;
const long double eps=1e-8;
double g[N][N];
double ans[N];
void gauss(){
	for(int i=1;i<=n;i++){
		int t=i;
		for(int j=i+1;j<=n;j++)
			if(g[t][i]<g[j][i])
				t=j;
		if(fabs(g[t][i])<eps){
			puts("No Solution");
			return ;
		}
		if(t!=i) swap(g[i],g[t]);
		
		double div=g[i][i];
		for(int j=i;j<=n+1;j++) g[i][j]/=div;
		for(int j=i+1;j<=n;j++){
			div=g[j][i];
			for(int k=i;k<=n+1;k++) g[j][k]=g[j][k]+(-div)*g[i][k];
		}
	}
	ans[n]=g[n][n+1];
	for(int i=n-1;i>=1;i--){
		for(int j=n;j>i;j--)
			g[i][n+1]-=ans[j]&g[i][j];
		ans[i]=g[i][n+1];
	}
	for(int i=1;i<=n;i++) printf("%.2lf\n",ans[i]);
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
			scanf("%lf",&g[i][j]);
	gauss();
	return 0;
}

这是判定解版的：

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n;
const int N=102;
const long double eps=1e-8;
double g[N][N];
int gauss(){
	int r,c;
	for(r=c=1;c<=n;c++){
		int t=r;
		for(int i=r+1;i<=n;i++)
			if(fabs(g[i][c])>fabs(g[t][c]))
				t=i;
		if(fabs(g[t][c])<eps) continue;
		if(t!=r) swap(g[t],g[r]);
		for(int i=n+1;i>=c;i--) g[r][i]/=g[r][c]; 
		for(int i=r+1;i<=n;i++)
			for(int j=n+1;j>=c;j--)
				g[i][j]-=g[r][j]*g[i][c];
		r++;
	}
	if(r<=n){
		for(int i=r;i<=n;i++) if(fabs(g[i][n+1])>eps) return 2;
		return 1;
	}
	for(int i=n-1;i>=1;i--)
		for(int j=i+1;j<=n;j++)
			g[i][n+1]-=g[j][n+1]*g[i][j];
	return 0;
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
		scanf("%lf",&g[i][j]);
	int op=gauss();
	if(op==0) for(int i=1;i<=n;i++) printf("%.2lf\n",g[i][n+1]);
	else if(op==1) puts("Infinite group solutions");
	else puts("No solution");
	return 0;
}

求解异或线性方程组

简述

异或线性方程组是长这样的：

\[\begin{cases}x_1\oplus x_2=1\\x_2\oplus x_3=0\\x_1=1\end{cases} \]

其增广矩阵：

\[\begin{bmatrix}1& 1& 0 &|&1\\0&1&1 &|&0\\1&0&0&|&1\end{bmatrix} \]

事实上，求解异或线性方程组和普通的线性方程组没有太大区别。只是解和系数都变成了 \(0\) 和 \(1\)，运算符号从 \(+\) 变成了 \(\oplus\)。多了一个圆圈

性质

异或其实有一个名字叫做不进位加法。它也有和四则运算一样的性质。

1. 两异或线性方程交换位置，解不变
2. 两异或线性方程相互异或，解不变。

同样的，我们想要构造出一个完美阶梯型矩阵。形如：

\[\begin{bmatrix}1& 0& 1 &|&1\\0&1&1 &|&0\\0&0&1&|&0\end{bmatrix} \]

过程

从左向右遍历每一列，我们找到一个该列系数为 \(1\) 的数，作为主方程。

然后将下面的方程该列为 \(1\) 的方程，和主方程相异或，就能把该列的系数全部化为 \(0\)。

判定解和回带操作也大同小异。

只有 g[i][j] 为 \(1\) 时，我们才能回带，则此时我们直接 g[i][n+1]^=g[i][j]&g[j][n+1] 就行了。

题目

高斯消元解异或线性方程组

点击查看代码

#include<bits/stdc++.h>
using namespace std;
int n;
const int N=102;
int g[N][N];
int gauss(){
	int r,c;
	for(r=c=1;c<=n;c++){
		int t=r;
		for(int i=r+1;i<=n;i++)
			if(g[i][c]){
				t=i;
				break;
			}
		if(!g[t][c]) continue;
		if(t!=r) swap(g[t],g[r]);
		for(int i=r+1;i<=n;i++)
			if(g[i][c])
				for(int j=c;j<=n+1;j++)
					g[i][j]^=g[r][j];
		r++; 
	}
	if(r<=n){
		for(int i=r;i<=n;i++) if(g[i][n+1]) return 2;
		return 1;
	}
	for(int i=n-1;i>=1;i--) 
		for(int j=i+1;j<=n;j++)
			g[i][n+1]^=g[i][j]*g[j][n+1];
	return 0;
}
int main(){
	scanf("%d",&n);
    for(int i=1;i<=n;i++)
		for(int j=1;j<=n+1;j++)
			scanf("%d",&g[i][j]);
	int op=gauss();
	if(op==0) for(int i=1;i<=n;i++) printf("%d\n",g[i][n+1]);
	else if(op==1) puts("Multiple sets of solutions");
	else puts("No solution");
	return 0;
}