入门到入土 | 高斯消元
在谈高斯消元之前，我们需要先了解一下高斯消元是干什么的，是求解什么的方法，知道矩阵和行列式的大佬可以直接跳到最后面
- PS：由于高斯消元是主要听HYL讲的，于是这篇博客也部分借鉴了HYL的原博客，可能我的描述方式和顺序会略有区别，有些见解可能更偏向于我这种初学者，因此部分言论可能并不严谨，甚至出现错误，请发现问题的大佬及时提出问题，Deltana在此感激不尽
- PPS：本篇博客也是我第一次用\(\LaTeX\),手生导致的格式爆炸警告！
高斯消元
简单来说，高斯消元是一种计算线性方程组的方法
那么说到线性方程组，就不得不提到行列式以及矩阵，而高斯消元正是建立在矩阵或行列式上的一种计算线性方程组的方法，那么我们一点一点来，先说矩阵和行列式
矩阵
什么是矩阵
在数学中，矩阵（Matrix）是一个按照长方阵列排列的复数或实数集合，最早来自于方程组的系数及常数所构成的方阵
比如说这样：
\(\begin{vmatrix} 1 & 2 \\ 3 & 4 \\8 & 9 \end{vmatrix}\)
这就是一个由实数元素构成的实矩阵
\(\begin{vmatrix}2i+1 & 2i+1 \\ -i & 7i\end{vmatrix}\)
这就是一个由复数元素构成的复矩阵
我们把第\(i\)行第\(j\)列的元素\(a_{ij}\)称为这个矩阵的\((i,j)\)元
我习惯把\(m*n\)的矩阵\(A\)表示为\(A_{mn}\)
下面来简单介绍一下矩阵的基本运算
矩阵的基本运算
由于我太菜了不会证明所以证明就先咕了
虽然绝大部分矩阵基本运算在下面的高斯消元中我们用不到，但是我还是写出来了（
线性运算
我们把以下介绍到的运算称为线性运算：
- 矩阵加减法
- 矩阵乘一个数
矩阵加减法
对于\(A_{mn}和B_{mn}\)，\(A+B=\)
\(\begin{bmatrix}a_{11} & a_{12} & a_{13} & \cdots & a_{1n}\\a_{21} & a_{22} & a_{23} & \cdots & a_{2n}\\\vdots & \vdots & \vdots & \ddots & \vdots\\a_{m1} & a_{m2} & a_{m3} & \cdots & a_{mn}\\\end{bmatrix} + \begin{bmatrix}b_{11} & b_{12} & b_{13} & \cdots & b_{1n}\\b_{21} & b_{22} & b_{23} & \cdots & b_{2n}\\\vdots & \vdots & \vdots & \ddots & \vdots\\b_{m1} & b_{m2} & b_{m3} & \cdots & b_{mn}\\\end{bmatrix} = \begin{bmatrix}a_{11}+b_{11} & a_{12}+b_{12} & a_{13}+b_{13} & \cdots & a_{1n}+b_{1n}\\a_{21}+b_{21} & a_{22}+b_{22} & a_{23}+b_{23} & \cdots & a_{2n}+b_{2n}\\\vdots & \vdots & \vdots & \ddots & \vdots\\a_{m1}+b_{m1} & a_{m2}+b_{m1} & a_{m3}+b_{m1} & \cdots & a_{mn}+b_{mn}\\\end{bmatrix}\)
相信也很明了了，也就是两个矩阵对应的\((i,j)\)元相加
矩阵加法满足：
- \(A+B=B+A\)
- \((A+B)+C=A+(B+C)\)
\(A-B\)同理，是对应的\((i,j)\)元相减
矩阵的数量乘法
请注意这并不是矩阵乘法，矩阵乘法下面会提到
对于\(\mu和A_{mn}\)，\(\mu*A_{mn}=\)
\(\mu*\begin{bmatrix}a_{11} & \cdots & a_{1n}\\\vdots & \ddots & \vdots\\ a_{m1} & \cdots & a_{mn}\end{bmatrix}=\begin{bmatrix}\mu*a_{11} & \cdots & \mu*a_{1n}\\\vdots & \ddots & \vdots\\ \mu*a_{m1} & \cdots & \mu*a_{mn}\end{bmatrix}\)
也就是让矩阵中的每一个元素都和\(\mu\)相乘
此计算满足：
-
\((\mu\nu)*A=\mu*(\nu A)\)
-
\((\mu+\nu)*A=\mu A+\nu A)\)
-
\((A+B)*\mu=\mu A+\mu B)\)
其它运算
转置
也就是将\(A_{m n}\)的第\(i\)行换成第\(i\)列
表示为\(A^T=\begin{bmatrix}a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \cdots & a_{mn}\end{bmatrix}^T=\begin{bmatrix}a_{11} & \cdots & a_{m1} \\ \vdots & \ddots & \vdots \\ a_{1n} & \cdots & a_{mn}\end{bmatrix}\)
转置运算满足：
- \((A ^ T) ^ T = A\)
- \((a A) ^ T = a A ^ T\)
- \((A B) ^ T = A^T B^T\)
共轭&共轭转置
这两个运算是针对存在复数元的矩阵的，由于我对复数了解较少就不多做解释
矩阵\(A\)的共轭矩阵表示为\(\overline{A}\)，共轭转置矩阵为\(A^*=\overline{(A)^T}=\overline{A^T}\)
设\(A=\begin{bmatrix}3+i & 6 \\ i & -i+5\end{bmatrix}\)
则\(\overline{A}=\begin{bmatrix}3-i & 6 \\ -i & i+5\end{bmatrix}\)
则\(\overline{A^T}=\begin{bmatrix}3+i & -i \\ 6 & i+5\end{bmatrix}\)
矩阵乘法运算
矩阵乘法仅有在相乘的第一个矩阵的列数等于第二个矩阵的行数时才有定义
也就是说：对于\(A_{mn}和B_{kp}\)，只有\(n=k\)时矩阵乘法才有定义
对于矩阵\(A_{mn}和B_{nk}\)相乘，结果矩阵为\(C_{mk}\)
\(C\)中任意一个元素\(c_{ij}=a_{i1} * b_{1i}+a_{i2}*b_{2i}+\cdots+a_{in}+b_{ni}=\sum\limits_{r=1}^n a_{ir}*b_{rj}\)
为了方便理解，我举个栗子：
\(\begin{bmatrix}2 & 0 \\ 7 & 7\end{bmatrix} * \begin{bmatrix}-1 & 0 & 2\\ 3 & 6 & 5\end{bmatrix}=\begin{bmatrix}-1*2 + 3*0 & 0*2 + 6*0 & 2*2 + 0*5 \\ -1*7 + 3*7 & 0*7 + 6*7 & 2*7 + 5*7 \end{bmatrix}\)
矩阵乘法运算满足：
- \((AB)C=A(BC)\)
- \((A+B)C=AC+BC\)
- \(C(A+B)=CA+BC\)
注意：\(AB≠BA\)
特殊矩阵
- 对角矩阵，也就是除了主对角线\((i=j)\)以外的元素/系数皆为0
- 三角矩阵，分为上三角和下三角，前者是左下方的元素/系数皆为0，后者是右上方
那么到这里矩阵的基础知识就差不多介绍完了，下面来看看行列式
行列式
N阶行列式由\(N^2\)个数\(a_{ij}(1≤i,j≤,n)\)组成
求和公式
\(\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}={\sum\limits_{j_1,j_2,\cdots,j_n}sgn(j1j2j3{\cdots}jn)*a_{1{j_1}}*a{2{j_2}}\cdots*a{n{j_n}}}\)
这里\(sgn\)函数为:如果该全排列为偶，则返回1，否则则返回-1
即\(sgn(j1j2j3⋯jn)=(−1)^{τ(j1,j2,...jn)}\)
其中j是一个1至n的全排列
也就是说，从每一行，每一列中都取出一个数来，相乘，把枚举所得到的所有结果相加，就是最终答案
行列式的性质
这里面有我们下面要用到的/kk
由于我太菜了不会证明所以证明我又双叒叕先咕了
这里给出HYL的这部分的链链接，她讲的挺详细的（）
[Ⅰ]行列式和它的转置行列式相等
这里的转置和我刚才在矩阵那里说的转置是一样的
也就是说\(\begin{vmatrix}a & b \\ c & d\end{vmatrix}=\begin{vmatrix}a & c \\ b & d\end{vmatrix}\)
[Ⅱ]对换行列式的两行（列），行列式变号
也就是说\(\begin{vmatrix}a & b \\ c & d\end{vmatrix} = -\begin{vmatrix}c & d \\ a & b\end{vmatrix} = -\begin{vmatrix}b & a \\ d & c\end{vmatrix}\)
[Ⅲ]若行列式有两行（列）完全相同或成比例，此行列式等于零
这个就不举例子了，简单易懂
[Ⅳ]用一个数乘行列式的某行（列）等于用此次数乘此行列式
即：\(k*\begin{vmatrix}a & c \\ b & d\end{vmatrix} = \begin{vmatrix}a & b \\ k*c & k*d\end{vmatrix}\)
推论：行列式的某行（列）的所有元素的公因数可以提出来，行列式值不变
即：\(\begin{vmatrix}a & b \\ k*c & k*d\end{vmatrix}=k*\begin{vmatrix}a & c \\ b & d\end{vmatrix}\)
这个有时候能起到简化运算的作用
（想起来被 手模7阶行列式值 支配的恐惧）
[Ⅴ]如果行列式中某一行等于两个行列式之和，则这个行列式等于两个行列式之和，这两个行列式分别以这两组数为改行，其余各行相同
即：
若\(a_{2i}=b_{2i}+c_{2i}\)，则有：
\({\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}={\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\b_{21}&b_{22}&\cdots&b_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}+{\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\c_{21}&c_{22}&\cdots&c_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}\)
[Ⅵ]把一行乘以某数加到另一行，行列式的值不变
这个就是高斯消元的精髓（）
即：
若有\(b_{2i}=a_{kj}*\mu\)，则有：
\({\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}={\begin{vmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\ b_{21}&b_{22}&\cdots&b_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{vmatrix}}\)
再谈高斯消元
嘤嘤嘤终于写到这了
前面说到过，高斯消元是我们求解线性方程组的一种方法，其实它本质上是我们计算行列式的一种方法，略做扩展后，它也可以在矩阵里消元，再利用矩阵来计算线性方程组
前面我有提到过，计算行列式是有公式的，但是这个东西的复杂度是\(O(n*n!)\)......很明显这个指数级时间直接炸裂开，那么我们怎么计算呢？
先引入一种特殊的行列式：三角行列式，它和三角矩阵的定义基本一致
\(\begin{vmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ 0 & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{nn} \end{vmatrix}\)
那么我们在求这个行列式的时候，如果你简单举个例子就会发现：这个行列式的值就是主对角线的乘积
再有我们上面说到的性质Ⅵ，实际上所有行列式都可以换成这种形式
也就是说,我们可以根据一定的操作之后，我们可以把求值的复杂度大大缩减
实际上高斯消元也就是完成的这个工作，我们通过最后一条性质就可以对整个行列式进行反复的消元，最终达到上三角的形式
代码实现1.0
让我们一步一步来，先考虑最简单的
即用主对角线的每一个元素来消去下面每一行对应的应消为0的元素，用\({z_{ji}} \over {b_{ii}}\)来计算每次前面的行要为了消元要乘上的数
于是乎我们就得到了这个最为朴素的算法的实现形式，代码中应该已经有了比较多的注释
Guass1.0
在放\(Code\)之前，为了缩短篇幅，我先打出我一边所用的头文件以及我的瞎\(define\)，在下面的代码中我将不再放出这个
#include <bits/stdc++.h>
#define Heriko return
#define Deltana 0
#define S signed
#define LL long long
#define R register
#define I inline
using namespace std;
I void fr(LL &x)
{
LL f=1;char c=getchar();
x=0;
while(c<'0'||c>'9')
{
if(c=='-') f=-1;
c=getchar();
}
while(c>='0'&&c<='9')
{
x=(x<<3)+(x<<1)+c-'0';
c=getchar();
}
x*=f;
}
实际上看到1.0就应该知道后面还没完了吧
double z[114][114];//二维数组模拟矩阵
int n;//n阶矩阵
I double Guass1()
{
double ans=1.00;//最终要输出的结果
for(R LL i=1;i<=n;i++)
{
for(R LL j=i+1;j<=n;j++)
{
double temp=z[j][i]/z[i][i];//计算要计算的比值
for(R LL u=1;u<=n;u++) z[j][u]-=temp*z[i][u];//把要消去的那一行用当前行消(减)
}
}
for(R LL i=1;i<=n;i++) ans*=z[i][i];//计算主对角线的乘积之和，也就是上三角时行列式的值
Heriko ans;//返回答案
}
根据这个代码我们可大约得出，这个算法的时间复杂度约为\(O(n^3)\)，虽然还是不小，但是相比于指数级还是小了很多
但是这个\(1.0\)的\(Code\)并不完美，或者说是有许多的漏洞，主要有以下两点：
- 若要消去的元素本来就是0时，这个算比值的方法就无法成立了
- 精度损失很大
那么下面对这份代码的这两个问题进行逐一修改
代码实现1.0优化
我们先来针对第一个问题进行修改
根据我们前面所讲到的性质Ⅱ“将行列式中的两行（列）对换，行列式变号”，我们可以把对应元素为不为0的那一行和当前行进行调换，然后将行列式取负即可
Guass1.0+
dd a[514][514];
int n;//同1.0
I double Guass1Plus()
{
double ans=1.00;
for(R LL i=1;i<=n;i++)
{
for(R LL j=i+1;j<=n;j++)
{
if(fabs(z[j][i])>1e-8)//如果当前行的对应元素不为0，把它交换上来
{
ans=-ans;//因为我们交换了两行，所以行列式变号
for(R LL u=1;u<=n;u++) swap(z[i][u],a[j][u]);//交换这两行
break;//当前行已经交换完了
}
}
if(fabs(z[i][i])<=1e-8) Heriko 0.00;//若还是0，说明这一列都是0，
//根据前面说到的性质Ⅲ，这个行列式值就为0
for(R LL j=i+1;j<=n;j++)
{
double temp=z[j][i]/z[i][i];
for(R LL u=1;u<=n;u++) z[j][u]-=temp*z[i][u];//和刚才1.0一样的进行消元
}
}
for(R LL i=1;i<=n;i++) ans*=z[i][i];//计算行列式的值
Heriko ans;
}
但是这样之后我们损失精度最严重的地方还是没有解决
下面我们来分析如何解决精度问题
再我们上面的程序当中，损失精度最多的来源于这样一句
double temp=z[j][i]/z[i][i];
那么我们就要着眼于这句话的优化，让我们来逐步地分析
Q1:为什么会有精度损失？
A1:作为一个由整数元素构成的行列式，其值也必定是整数，但是我们这样算必定得出的时实数，而实数在运算的过程当中就不可避免的会出现精度损失
Q2:为什么这个精度损失大？
A2:因为我们没有做任何的处理，直接就让两个元素求比值，这样做的\(z[i][i]\)的大小是不确定的，精度损失也就相应的很大
Q3:如何解决？
A3:在保持这种算法基本不变的情况下，我们只能将损失降到尽可能小，因为对于double/longdouble而言，它在内存中的占用空间是一定的，也就是说它所存的数的位数是一定的，假如我们把所有的位数都给小数位去用，那么自然精度就会高的多。因此，为了使更多的空间被使用在小数部分上，我们要让这个比值尽可能的小，即让\({z[j][i]} \over {z[i][i]}\)尽可能的小，那么我们就需要让它的分母尽可能的大，也就是说让\(z[i][i]\)尽可能大
因此我们若想在不改变算法的前提下优化精度问题，就想办法把每列绝对值大的弄到上面来
Guass1.0++
据说这个貌似就能过掉洛谷的模板题了
I double Guass1PlusPlus()
{
double ans=1.00;
for(R LL i=1;i<=n;i++)
{
for(R LL j=i+1;j<=n;j++)
{
if(fabs(z[j][i])>fabs(z[i][i]))//若当前j遍历到的这个元素绝对值较大，交换两行
{
ans=-ans;//交换两行行列式变号
for(R LL u=1;u<=n;u++) swap(z[a][c],z[b][c]);
}
if(fabs(z[i][i])<=1e-8) Heriko 0.00;//若还是0，说明此列都为0，那么此行列式为0
for(R LL j=i+1;j<=n;j++)
{
double temp=z[j][i]/z[i][i];
for(R LL u=1;u<=n;u++) z[j][i]-=temp*z[i][i];
}
}
}
for(R LL i=1;i<=n;i++) ans*=z[i][i];
Heriko ans;
}
那么这个解决了两个问题后的这种算法叫做主元高斯消元法
但是由于算法的先天性劣势，我们无法完全消除精度损失
看来如果我们不改变实现算法是无法消除精度损失的了，为了解决掉这个问题，下面将引出辗转相消高斯消元法
代码实现2.0
那么我们先要理解这个算法，从这个名字入手
辗转相消，这个仿佛在哪听过......
对了！是辗转相除求最大公因数！那我们先来看看我们常用的GCD是怎么求的吧！
I LL GCD(R LL x, R LL y)
{
Heriko x==0 ? y : GCD(y,x%y);
}
我们也要参考辗转相除的方法来辗转相消，把下面那个三角全部消成0！
Guass2.0
I int Guass2()
{
int ans=1;
for(R int i=1;i<=n;i++)
{
for(R int j=1;j<=n;j++)
{
while(z[j][i]!=0)
{
int temp=z[i][i]/z[j][i];
for(R int u=1;u<=n;u++) z[i][u]-=z[j][u]*temp;
for(R int u=1;u<=n;u++) swap(z[i][u],z[j][u]);
ans=-ans;
}
}
}
for(R int i=1;i<=n;i++) ans*=z[i][i];
Heriko ans;
}
由于我不是很会推这个的时间复杂度，这里说一下结果就好了，是\(O(n^3)\)
具体的推导可以看HYL的推导过程
由此我们就实现了高斯消元，但是这只是建立在行列式上面的，如何用它来解决线性方程组呢？
利用高斯消元求线性方程组
刚才说到的洛谷模板就是求这个,题面推荐看前方灰色链接
刚才我们提到的是对于\(n*n\)行列式的高斯消元，实际上对于一个\(n*n+1\)矩阵应该也是适用的，于是乎我们就可以把非0次项系数和0次项系数（常数）都扔进这个n*n+1的矩阵！
于是乎我们要消元的就是这个矩阵
\(\begin{bmatrix}a_{11} & \cdots & a_{1n+1} \\ \vdots & \ddots & \vdots \\ a_{n1} & \cdots & a_{nn+1}\end{bmatrix}\)
消时候，我们不用管最右面那比行列式多出来的一列，按照\(n*n\)的矩阵消元就可，消完之后我们就可以发现实际上得到的是这样的一个矩阵
\({\begin{bmatrix}a_{11}&a_{12}&\cdots&a_{1n}&a_{1n+1}\\0&a_{22}&\cdots&a_{2n}&a_{2n+1}\\\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&\cdots&a_{nn}&a_{nn+1}\end{bmatrix}}\)
因为我们扔进去的元素是原方程组的系数，所以实际上最后一行变为了一个一元一次方程，这个方程的\(x_n\)我们是可以求解的，如何把它当作已知条件带回，那么每个行那个方程都变成一元一次方程，我们就可以求解最终的\(x\)
P3389 Code
#include <bits/stdc++.h>
#define Heriko return
#define Deltana 0
#define S signed
#define R register
#define I inline
#define D double
#define ll long long
#define ull unsigned long long
#define N 102
#define M number
using namespace std;
int n;
D z[N][N],ans[N];
inline bool guass(){
for(int i=1;i<=n;i++){
for(int j=i+1;j<=n;j++)
if(fabs(z[j][i])>fabs(z[i][i]))
for(int u=1;u<=n+1;u++) swap(z[i][u],z[j][u]);
if(fabs(z[i][i])<=1e-8) Heriko Deltana;
for(int j=i+1;j<=n;j++){
D k=z[j][i]/z[i][i];
for(int u=1;u<=n+1;u++) z[j][u]-=z[i][u]*k;
}
}
Heriko 1;
}
S main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)for(int j=1;j<=n+1;j++)scanf("%lf",&z[i][j]);
if(!guass()){printf("No Solution");Heriko 0;}
for(int i=n;i>=1;i--)
{
D x=0.00;
for(int j=i+1;j<=n;j++){
x+=ans[j]*z[i][j];
}
ans[i]=(z[i][n+1]-x)/z[i][i];
}
for(int i=1;i<=n;i++){
printf("%0.2lf\n",ans[i]);
}
Heriko Deltana;
}
End
我可算是写完了，整整搞了一天
都过二十四点了，好困，睡觉去了（
