高斯消元

前言

不会线性代数。

线性方程组

一个线性方程组通常是如下的形式：

\[\large \left\{\begin{matrix} 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_{2n}x_n = b_2\\ \vdots\\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n = b_m\\ \end{matrix}\right. \]

为了简便，下面会表示为：

\[\large \left\{\begin{matrix} &a_{1,1} &a_{1,2} &\cdots &a_{1,n - 1} &a_{1,n} | &a_{1,n + 1}\\ &a_{2,1} &a_{2,2} &\cdots &a_{2,n - 1} &a_{2,n} | &a_{2,n + 1}\\ &&&&\vdots\\ &a_{m,1} &a_{m,2} &\cdots &a_{m,n - 1} &a_{m,n} | &a_{m,n + 1}\\ \end{matrix}\right. \]

最简单的消元法自然是加减消元。
比如现在要把一个未知数 \(x_1\) 消去，那么就取出方程组里的一个方程 :

\[\large a_1x_1 + a_2x_2 + \cdots + a_nx_n = p_1 \]

方程的两侧同除 \(a_1\) 得到：

\[\large x_1 + \frac{a_2}{a_1}x_2 + \cdots + \frac{a_n}{a_1}x_n = \frac{p_1}{a_1} \tag{1} \]

现在就可以把这个方程代入其他方程，例如我们再找一个：

\[\large b_1x_1 + b_2x_2 + \cdots + b_nx_n = p_2 \tag{2} \]

把 \((1)\) 扩大 \(b_1\) 倍，用 \((2)\) 减去 \((1)\) 就消去了 \((2)\) 中的 \(x_1\)，这就是小学生都会的最简单的加减消元的过程。

然后可以发现在得到式子 \((1)\) 的时候有一个同除 \(a_1\) 的操作，如果这个 \(a_1\) 太小就会导致浮点数精度误差，于是我们每次从剩下的方程中选择主元系数（这里需要的是绝对值最大）最大的方程组来减小误差。

可以发现这个方程最后的形式会从：

\[\large \left\{\begin{matrix} &a_{1,1} &a_{1,2} &\cdots &a_{1,n - 1} &a_{1,n} | &a_{1,n + 1}\\ &a_{2,1} &a_{2,2} &\cdots &a_{2,n - 1} &a_{2,n} | &a_{2,n + 1}\\ &&&&\vdots\\ &a_{m,1} &a_{m,2} &\cdots &a_{m,n - 1} &a_{m,n} | &a_{m,n + 1}\\ \end{matrix}\right. \]

变成：

\[\large \left\{\begin{matrix} &a'_{1,1} &a'_{1,2} &\cdots &a'_{1,n - 1} &a'_{1,n} | &a'_{1,n + 1}\\ &0 &a'_{2,2} &\cdots &a'_{2,n - 1} &a'_{2,n} | &a'_{2,n + 1}\\ &&&&\vdots\\ &0 &0 &\cdots &0 &a'_{m,n} | &a'_{m,n + 1}\\ \end{matrix}\right. \]

然后就是倒推一下，\(x_n = \dfrac{a'_{m,n + 1}}{a'_{m,n}}\)，再从下往上把 \(x_n,x_{n - 1},\dots\) 都求出来即可。

如果某一行前 \(n\) 个数为 \(0\) 但是第 \(n + 1\) 个数不为 \(0\)，那么方程组无解。
如果某一行每个数都是 \(0\) 则有无穷多解。

int GaussElimin(int n,int m,db *res) {
	int r = 1,c = 1;
	for(;r <= n && c <= m;++c,++r) {
		int mx = r;
		rep(i,r + 1,n) if(std::abs(eq[i][c]) > std::abs(eq[mx][c]))
			mx = i;
		std::swap(eq[r],eq[mx]);
		if(std::abs(eq[r][c]) < eps) {
			--r;
			continue;
		}
		rep(i,r + 1,n) if(std::abs(eq[i][c]) > eps) {
			db t = eq[i][c] / eq[r][c];
			rep(j,c,m + 1)
				eq[i][j] -= eq[r][j] * t;
			eq[i][c] = 0;
		}
	}
	rep(i,r,m) if(std::abs(eq[i][c]) > eps)
		return -1;// No Solution
	if(r <= m) return m - r + 1;//Free Variables Exist
	repb(i,n,1) {
		rep(j,i + 1,m)
			eq[i][m + 1] -= eq[i][j] * res[j];
		res[i] = eq[i][m + 1] / eq[i][i];
	}
	return 0;
}

然后复杂度也很显然，是 \(\mathcal{O}(n^3)\) 的。

这时候天降了一个名叫 约旦 (Jordan，又译若尔当)的猛男，他说这个消元还需要回代，有一些 \(naive\)。如果能把矩阵最后转化成：

\[\large \left\{\begin{matrix} &1 &0 &\cdots &0 &0 | &a'_{1,n + 1}\\ &0 &1 &\cdots &0 &0 | &a'_{2,n + 1}\\ &&&&\vdots\\ &0 &0 &\cdots &0 &1 | &a'_{m,n + 1}\\ \end{matrix}\right. \]

那就直接出答案了。

于是他发明了 高斯-约旦消元，缺点是不能判断一个“无解”是方程无解还是方程有无数解（有自由元）。

int GaussJordan(int n,db *res) {
	rep(i,1,n) {
		int mx = i;
		rep(j,i + 1,n) if(std::abs(eq[j][i]) > std::abs(eq[mx][i]))
			mx = j;
		if(std::abs(eq[mx][i]) < eps)
			return 0;
		std::swap(eq[i],eq[mx]);
		db t = eq[i][i];
		rep(j,1,n + 1)
			eq[i][j] /= t;
		rep(j,1,n) if(i != j) {
			t = eq[j][i];
			rep(k,1,n + 1)
				eq[j][k] -= t * eq[i][k];
		}
	}
	rep(i,1,n)
		res[i] = eq[i][n + 1];
	return 1;
}

高斯消元解线性方程组应用

球形空间产生器

给出 \(n\) 维空间一个球上 \(n + 1\) 个点坐标，求球心。
首先可知每个点到球心距离相等，令球心坐标为 \((x_1,x_2,x_3,\cdots,x_n)\)，第 \(i\) 个点第 \(j\) 维的坐标为 \(p_{i,j}\)，那么有 :

\[\large \forall i \in[1,n+1],\sum_{j = 1}^{n} (p_{i,j} - x_j)^2 = r^2 \]

于是现在有 \(x_1\sim x_n\) 和 \(r\) 一共 \(n + 1\) 个未知量以及 \(n + 1\) 个等式，但是不是线性方程组，但是没有关系，把所有的式子排列起来然后：

\[\large \]

按照这个套路就消掉 \(r^2\) 了。
然后继续尝试摘掉平方，移项一下。

咕咕咕！

矩阵

现在我们从线性代数的角度理解高斯消元。

逆矩阵

矩阵行列式

异或方程组