高斯消元算法
高斯消元其实在算法竞赛中算是一个十分常见的算法。它的大致思想就和初中阶段学到的加减消元法差不多。这个算法的时间复杂度为\(O(n^3)\)，是一个相当简单的算法，但是具体实现需要一些思考。
1.1 问题引入
给定方程组\(\begin{cases}x+3y+4z=5\quad(1)\\x+4y+7z=3\quad(2)\\9x+3y+2z=2\quad(3)\end{cases}\),求解之。
这是一个相当简单的三元一次方程组，直接用加减消元就可以得出解。当然，这里给出一个比较不讨巧的消元顺序，为了方便后面的理解。
首先（2）式减去（1）式：\(y+3z = -2\quad(2)'\)
接下来\(9\times(1)-(3)\):\(24y+34z=43\quad(3)'\)
此时\(y\)和\(z\)就构成了二元一次方程组了。
计算\(24\times(2)'-(3)'\): \(38z=-91\).
解得\(z=-\frac{91}{38}\).此时不断的回代可以得到\(y=\frac{197}{38}\), \(x=-\frac{37}{38}\).
我们可以把整个计算过程得到的新的方程式列出来：
仔细观察一下这个方程组，这对于之后的解题有帮助。
思考一下，对于任意的三元一次方程组，或者更进一步的，对于任意的\(n\)元一次方程组，我们能不能设计一个通用的算法，求解方程组呢？
1.2 初等列变换
我们把\(x\),\(y\),\(z\)的系数提取出来，将它们列成一个3*3的表：\(\begin{bmatrix}1&3&4\\1&4&7\\9&3&2\end{bmatrix}\)
我们还可以把方程的右侧常数列在右边：\(\begin{bmatrix}1&3&4&|5\\1&4&7&|3\\9&3&2&|2\end{bmatrix}\)
这个3*4的表格叫做方程组的增广矩阵。规定矩阵的第\(i\)行为\(r_i\),如\(r_1=\begin{bmatrix}1&3&4&|5\end{bmatrix}\).
我们可以对它进行这样几个操作：
- 用一个非零常数乘上某一行。这个操作记为\(r_i=cr_i\)，其中\(c\)为非零常数
- 把其中一行的若干倍加至另一行上。这个操作记为\(r_i=cr_i+c'r_j\)
- 交换两行的位置。直接记作\(\text{swap}(r_i,r_j)\)就好了。
这样一来，我们之前的操作可以这样表示：
1.\(r_2=r_2-r_1\)
2.\(r_3=9r_1-r_3\)
3.\(r_3=24r_2-r_3\)
这个时候增广矩阵变成了这个样子：\(\begin{bmatrix}1&3&4&|5\\0&1&3&|-2\\0&0&38&|-91\end{bmatrix}\)
往往矩阵元素等于0时我们省略不写，上面的矩阵还可以这样写：\(\begin{bmatrix}1&3&4&|5\\ &1&3&|-2\\ & &38&|-91\end{bmatrix}\)
排版有点丑，但是我们还是可以看出来，这个矩阵左边的形状有点像一个倒过来的阶梯。人们就叫它简化阶梯型矩阵。把原增广矩阵变为现在的简化阶梯型矩阵的过程就叫做高斯消元。
算出了简化阶梯型矩阵之后，直接一步一步回代就可以了。
1.3 基本实现
总结一下，高斯消元的思路：
对于每一行\(r_i\)，我们要让\(r_{i,i}=1\)，并且\(r_{i,i}\)下面的元素均变为\(0\)。我们这样做：
- 选取\(|r_{j,i}|\)最大的一行\(r_j\)，将它与\(r_i\)交换。之后\(r_{i,i}\)就会变为当前列最大的元素。
- 把\(r_{i,i}\)变为\(1\)。很简单，直接让\(r_i=\frac{1}{r_{i,i}}r_i\)就行了。
- 把\(r_{i,i}\)下面的元素都变为零。对于\(i+1 \leq j \leq n\)，我们让\(r_j=r_j-\frac{r_{j,i}}{r_{i,i}}r_i\)，调整\(r_i\)的比例并把\(r_{j,i}\)消去。因为\(r_{i,i}=1\)，这个操作可以写成\(r_j=r_j-r_{j,i}r_i\)就行了。
这一段的大概代码如下：
RP(i,1,n)
{
rg int cur=i;
RP(j,1,n)
if(fabs(r[cur][i])<fabs(r[j][i]))
cur=j;
if(i!=cur)
swap(r[cur],r[i]);
db div=r[i][i];
RP(j,i,n+1)
r[i][j]/=div;
RP(row,i+1,n)
{
db rat=r[row][i];
RP(col,i,n+1)
r[row][col]-=rat*r[i][col];
}
}
在最后回代的过程中，由于最后一行的形式一定是一个一元一次方程，最后一行的解就是\(x_n=r_{n+1}\)。
之后的每一行\(r_i\):\(\sum\limits_{j=i}^n{x_j\cdot r_{i,j}}=r_{i,n+1}\)，移项就得到了解：$$x_i=\dfrac{r_{i,n+1}-\sum\limits_{j=i+1}^n{x_j\cdot r_{i,j}}}{r_{i,i}}=r_{i,n+1}-\sum\limits_{j=i+1}^n{x_j\cdot r_{i,j}}$$
x[n]=r[n][n+1];
DRP(i,n-1,1)
{
x[i]=r[i][n+1];
RP(j,i+1,n)
x[i]-=r[i][j]*x[j];
}
1.4 一些特判细节和其他注意事项
如果方程的某一列元素均为\(0\)，而常数项\(\neq 0\)，那么方程组是一定无解的。这里直接特判就可以了。
另外，在交换行时可以特判一下两行是否相等。减掉没必要的操作可以稍微提速。
这个算法还有一个加强版本:高斯-约旦消元法。这个算法其实就是把回代的操作直接用初等行变换实现，不过这个算法在矩阵求逆上有着很好的应用。
