【学习笔记】高斯消元法

HyperSQ · 2022-02-19 14:49:53 · 个人记录

高斯消元

一、线性方程组

求解形如

的 n 元一次方程组（也就是线性方程组），我们会用到高斯消元法。

我们把一个线性方程组简便记成上面这样的形式。

高斯消元主要依据线性方程组的以下三个性质：

1、方程组中的两条等式交换位置，解不变。

2、方程组中某个等式乘上某个非零实数，解不变。

3、方程组中某个等式减去另一个等式，解不变。

接下来我们通过一个实例来讲解高斯消元法。下文用 r 表示行。

由于一个纵列对应的是一个元，我们每次从左到右处理每个元。

对于第一个元，我们找到系数绝对值最大的那一项，并将其交换至第一列。（下文会解释这个步骤）

将第一列的第一个元系数化为 1 。同时，将其他行用性质3减去第一行的 k 倍，将其它式子中的第一个系数全部化为 0 。

发现右下角有一个子方阵，于是我们仿照上面的方法继续处理第二元，

最后将第三个元系数化为 1 ，

这样子，我们将原方程组化为了上三角形式。

化为上三角形式方便我们对解进行回代，接下来我们进行回代操作。

首先可以得到 x_3=-33 。我们把这个结果代入 r_2 可以得到 x_2=14。

再将 x_3 与 x_2 代入 r_1，我们又能解得 x_1=26。

于是方程组的解为

通过上面这个例子，聪明的你应该能够理解高斯消元算法的步骤了：

一、从 1 到 n 依次消去第 i 个元

二、消第 i 个元时，选择第 i 列中系数绝对值最大的未处理行 r，交换到第 i 行。

三、将第 i 行的主元系数化为 1。

四、将未处理行主元系数化为 0 。

五、将第 i 行看作是已处理行。继续对第 i+1 列进行消元。

通过上面的步骤，我们就能将原方程化为上三角的形式。再对元进行回代，就能得出方程组的解。

在步骤二中，如果行 r 主元系数依然为 0 ，则说明该方程组不存在唯一解，因为该元可取任意值。

如果最后的上三角形式中，存在一行系数全部是 0 而常数项不为 0 ，则方程组 无解。

一般来说，我们只需要判断是否存在唯一解即可。

这个算法时间复杂度为 O(n^3)，详见代码。

#include <bits/stdc++.h>
using namespace std;

const int maxn=2e2+5;
const double eps=1e-10;
int n;
double a[maxn][maxn];
double x[maxn];

int main(){
    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 r=i;
        for(int j=i+1;j<=n;j++){
            if(abs(a[j][i])>abs(a[r][i])){
                r=j;
            }
        }//寻找最大系数
        if(abs(a[r][i])<eps){//判断是否有解
            puts("No Solution");
            return 0;
        }
        if(r!=i)
        for(int j=1;j<=n+1;j++){
            swap(a[r][j],a[i][j]);
        }//交换
        double div=a[i][i];
        for(int j=i;j<=n+1;j++){
            a[i][j]/=div;
        }//系数化为一
        for(int j=i+1;j<=n;j++){
            double mul=a[j][i];
            for(int k=i;k<=n+1;k++){
                a[j][k]-=a[i][k]*mul;
            }
        }//消元
    }//通过这段程序，已经化为上三角形式
    x[n]=a[n][n+1];
    for(int i=n-1;i>=1;i--){
        x[i]=a[i][n+1];
        for(int j=i+1;j<=n;j++){
            x[i]-=a[i][j]*x[j];
        }
    }//回代
    for(int i=1;i<=n;i++){
        printf("%.2f\n",x[i]);
    }
}

简单说一说为什么要找绝对值最大的那一行，一是因为可以寻找一个非零行，二是作浮点数除法时，显然除以一个大数会让精度更好，误差更小。

二、高斯消元求行列式

由于行列式有特殊性质，所以解方程组的方法对行列式求解也适用。

我们知道，当行列式化为上三角形式时，行列式的值就是对角线的累积。

注意到我们在消元时作了整行除以 k 的操作，根据行列式的性质，相当于行列式的值也除以了 k ，于是我们要把 k 乘到最终答案里弥补损失。

另外，如果交换行时，根据行列式的性质，答案应该进行取反。

三、矩阵求逆

与线性方程组类似，矩阵可以进行交换、某行乘 k、某行加上另一行的操作。这称为矩阵的初等行变换。

实际上，这三个操作本质上分别对应着一种线性变换。

同样以上面那个矩阵为例，

交换第一第二行对应的线性变换是

第二行乘上 \frac{1}{2} 对应的线性变换是

第一行加上第二行的三倍对应的线性变换是

综上所述，我们对矩阵进行的三类操作其实均可看作是一次线性变换。

由逆矩阵的定义可知，

记 P=A^{-1}，说明如果可逆，那么 P 相当于就是对 A 进行了若干次初等行变换，并让 A 转变为 E 。即

对 A 的这些行变换本质上就是高斯消元的过程。

其中 P_n\cdots P_3P_2P_1 就等于 P，也就是 A 的逆矩阵。稍微作一下变换，

也就是说，在我们将 A 高斯消元的同时，对 E 也作同样的行变换，最后得到的“记录”就是 A 的逆矩阵。

套用高斯消元的模板即可。这里的逆矩阵是取模意义下的，所以所有的除法全部用乘法逆元代替。

具体细节就看代码吧。

#include <bits/stdc++.h>
#define ll long long
using namespace std;

const ll mod=1e9+7;
const int maxn=4e2+5;
ll a[maxn][maxn];
ll b[maxn][maxn];
int n;

ll qpow(ll a,ll b){
    ll ret=1;
    while(b){
        if(b&1) ret=ret*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ret;
}

inline ll inv(ll x){
    return qpow(x,mod-2);
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%lld",&a[i][j]);
    for(int i=1;i<=n;i++) b[i][i]=1;//初始化一个b矩阵为E，记录A中的操作
    //接下来要注意，B的记录操作一定要是整行变换的。
    //在一般的高斯消元中，A不进行整行操作，是因为行前已经削为0，但本质上还是进行了整行变换
    //但是B中行的数字都不一定是非0的，如果不进行整行变换就会出错。
    for(int i=1;i<=n;i++){
        int r=i;
        for(int j=i;j<=n;j++){//这里找一个非零值就可以了，因为整数运算不存在误差
            if(a[j][i]){
                r=j;break;
            }
        }
        if(!a[r][i]){//对角线中有0，即矩阵行列式为0，故不存在逆矩阵
            puts("No Solution");return 0;
        }
        if(r!=i) swap(a[r],a[i]),swap(b[r],b[i]);
        ll div=inv(a[i][i]);
        for(int j=1;j<=n;j++){
            a[i][j]=a[i][j]*div%mod;b[i][j]=b[i][j]*div%mod;
        }//注意B中整行的操作
        for(int j=i+1;j<=n;j++){
            ll mul=a[j][i];
            for(int k=1;k<=n;k++){//还是要整行
                a[j][k]-=mul*a[i][k]%mod;
                a[j][k]=(a[j][k]+mod)%mod;
                b[j][k]-=mul*b[i][k]%mod;
                b[j][k]=(b[j][k]+mod)%mod;
            }
        }
    }
    for(int i=n-1;i>=1;i--){//最后的回代过程还是要注意整行相减！！整行！！所以是O(n^3)的
        for(int j=i+1;j<=n;j++){
            for(int k=1;k<=n;k++){
                b[i][k]-=b[j][k]*a[i][j]%mod;
                b[i][k]=(b[i][k]+mod)%mod;
            }
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            printf("%lld ",b[i][j]);
        }
        puts("");
    }
}

高斯消元还有很多应用，日后掌握了再补充。