浅谈高斯消元法

永远爱循环 · 2025-02-08 11:51:26 · 算法·理论

高斯消元法

1.介绍

高斯消元，用于求解m个n元一次方程组成的方程组。

2.操作

我们可以将一个线性方程组写成一个 m 行 n+1 列的增广矩阵。

同时，我们定义矩阵的初等行变换：

• 用一个非零的数乘以某一行
• 把其中一行的若干倍加到零一行上
• 交换两行的位置

我们可以对增广矩阵进行若干次初等行变换来求解方程组：

实际上矩阵还可以继续化简：

像这样通过初等行变换，把增广矩阵变为阶梯型矩阵就是高斯消 元算法。该算法的思想是，对于每一个未知量 x_i，找到一个 x_i 系数不 为零但是 x_1 ... x_{i−1} 系数都为零的方程，利用该方程通过初等行变换把其 它方程的 x_i 系数全部消成 0。

高斯消元

其中，若把增广矩阵转化为行梯阵式，如

\begin{bmatrix} 1 &0 &0 &1\ 0 &1 &0& -2\ 0 &0& 1 &3\ \end{bmatrix}

\begin{bmatrix} 1 &2 &−1 &−6\ 0 &1 &1& 1\ 0 &0& 1 &3\ \end{bmatrix}的意义是\begin{cases} &x_1+2x_2+(-1)x_3=-6\ &x_2+x_3=1\ &x_3=3\ \end{cases}，最后我们需要将x_3=6代回x_2+x_3=6这个方程，从而求解出x_2=-2，再把\begin{cases} &x_2=3\ &x_3=-2\ \end{cases}

3. 有解、无解、无数个解

无解

在高斯消元的过程中，有可能会出现 0=d 这样的方程，其中 d 是一 个非零常量。出现这种情况则意味着方程无解。

无数个解

其次，有可能找不到 x_i 系数不为零的方程。例如：

在这个例子中找不到 x_2 系数为 0 的方程，解可以写为

不论x_2取任何值，都可以求出对应的 x_1，也就是方程组有无穷多的解。这种情况下，我们把 x_1, x_3 称为主元，把 x_2 称为自由元。

判断

可以发现，在最后的阶梯型系数矩阵中（即高斯-约旦消元后的简化行梯阵式），每个主元仅仅在一个位置上 的系数 (i,j) 非零，且第 i 行、第 j 列其他位置都为 0。

• 若系数a_{i,j}为零：
  • 若b_i为0，可以判断方程组无解。
  • 否则，b_i不为0，可以判断方程组存在无穷多的解，即无数个解。
• 否则，系数a_{i,j}不为零，并且可知第i个线性方程其他系数为0，所以第i个未知数x_i=\frac{b_i}{a_{i,j}}，即有解。

注意，在判断此处的时候，一定是先判断是否无解，再判断是否有无数个解，例如：

我们如果先判断存在自由元为无数个解，那么我们极可能容易忽略接下来的线性方程自由元是否存在，即无解的情况。

4.代码实现

时间复杂度为O(N^3)。

• 高斯消元

  void Gauss(){
  for(int i=1;i<=n;i++){
      int pos=i;
      for(int j=i+1;j<=n;j++){
          if(fabs(a[j][i])>fabs(a[pos][i]))pos=j;
      }
      //在 i~n 之间找到一个绝对值较大的一行
      for(int j=1;j<=n+1;j++)swap(a[i][j],a[pos][j]);
      //交换
      if(fabs(a[i][i])<eps){
          //这里可以用 aii 代表 i~n 行的 xi 的系数是因为已经把绝对值最大的换到前面来了
          //如果绝对值最大的都为0，说明全是0
          //由于是 double 浮点数，所以会存在精度问题，所以用 eps=1e-6 代替0
  
          //因为第 i 行的主元为 xi ，如果它的系数为0,就变成了自由元
          //所以要么无解，要么无数个解
          puts("No Solution");
          return ;
      }
      double pre=a[i][i];
      //一定要先存下来，因为在 xi 的系数 aii 会变为1，之后其他的数就除以更改过的 aii 而不是原先的 aii
      for(int j=1;j<=n+1;j++) a[i][j]=1.0*a[i][j]/pre; //把 xi 的系数变为1
      for(int j=1+i;j<=n;j++){
          double rate=a[j][i];
          //同理要先存下来
          for(int k=1;k<=n+1;k++){
              a[j][k]=a[j][k]-rate*a[i][k];
              //将其他的 xi 的系数变为0
          }
      }
  }
  ans[n]=a[n][n+1];
  for(int i=n-1;i>=1;i--){
      ans[i]=a[i][n+1];
      for(int j=i+1;j<=n;j++)//之前的都为0，故不需要在处理
          ans[i]-=(a[i][j]*ans[j]);
  }//回带操作
  for(int i=1;i<=n;i++){
      printf("%.2lf\n",ans[i]);
      //保留两位小数输出
  }
  }

• 高斯-约旦消元

  void Gauss(){
  for(int i=1;i<=n;i++){
      int pos=i;
      for(int j=i+1;j<=n;j++){
          if(fabs(a[j][i])>fabs(a[pos][i]))pos=j;
      }
      //在 i~n 之间找到一个绝对值较大的一行
      for(int j=1;j<=n+1;j++)swap(a[i][j],a[pos][j]);
      //交换
      if(fabs(a[i][i])<eps){
          //这里可以用 aii 代表 i~n 行的 xi 的系数是因为已经把绝对值最大的换到前面来了
          //如果绝对值最大的都为0，说明全是0
          //由于是 double 浮点数，所以会存在精度问题，所以用 eps=1e-6 代替0
  
          //因为第 i 行的主元为 xi ，如果它的系数为0,就变成了自由元
          //所以要么无解，要么无数个解
          puts("No Solution");
          return ;
      }
      double pre=a[i][i];
      //一定要先存下来，因为在 xi 的系数 aii 会变为1，之后其他的数就除以更改过的 aii 而不是原先的 aii
      for(int j=1;j<=n+1;j++) a[i][j]=1.0*a[i][j]/pre; //把 xi 的系数变为1
      for(int j=1;j<=n;j++){
          if(i==j) continue;
          double rate=a[j][i];
          //同理要先存下来
          for(int k=1;k<=n+1;k++){
              a[j][k]=a[j][k]-rate*a[i][k];
              //将其他的 xi 的系数变为0
          }
      }
  }
  for(int i=1;i<=n;i++){
      printf("%.2lf\n",a[i][n+1]);
      //保留两位小数输出
  }
  }

5.高斯消元法的拓展使用

• 逆矩阵

  • 介绍

  定义一个n*n的单位矩阵:

  I_n= \begin{bmatrix} 1 & 0& 0 &5\\ 0 & 0 &0 &0\\ 0 &0 & 0 &2\\ \end{bmatrix}

  对于一个n*m的矩阵A来说，若存在一个m*n的矩阵B，满足AB=I_n，BA=I_m，则称B为A的逆矩阵，记作A^{-1}。实际上，如果矩阵A可逆，那么A^{-1}唯一，且n=m。

  • 性质

  $3.$转置的逆：如果矩阵$A$可逆，那么矩阵$A$的转置矩阵$A^T$也是可逆的，并且$(A^T)^{-1}=(A^{-1})^T$。 $4.$逆矩阵的逆：如果矩阵$A$可逆，那么矩阵$A$的逆矩阵的逆矩阵为矩阵$A$，即$(A^{-1})^{-1}=A$。 + #### 求法 注：此处谈及逆矩阵的目的是为了拓展使用高斯消元，所以其他方法并不详细介绍。 $1.$待定系数法。 $2.$伴随矩阵求逆矩阵 $3.$初等变换求逆矩阵： $\text{对于矩阵} A = \begin{bmatrix} 2& −1 &0\\ −1& 2 &−1\\ 0 &−1 &2\\ \end{bmatrix} ，\text{我们把它和 }I_3 \text{拼在一起，得到:}\\ [A|I]= \left[\begin{array}{ccc|ccc} 2 & -1 & 0 & 1 & 0 & 0 \\ -1 & 2 & -1 & 0 & 1 & 0 \\ 0 & -1 & 2 & 0 & 0 & 1 \end{array}\right]\\ \text{然后通过初等行变换把}A \text{变为单位矩阵：}\\ [I|B]=\left[\begin{array}{ccc|ccc} 1 & 0 & 0 & 0.75 & 0.5 & 0.25 \\ 0 & 1 & 0 & 0.5 & 1 & 0.5 \\ 0 & 0 & 1 & 0.25 & 0.5 & 0.75 \end{array}\right]\\ \text{右边的}B \text{矩阵即为 }A^{−1}\text{。}

  • 代码实现

    注意：建议使用高斯-约旦消元，高斯消元回代时实现代码复杂。

    
    void Gauss(){//高斯消元
    for(int i=1;i<=n;++i){
      int pos=i;
      for(int j=i+1;j<=n;++j){
          if(a[j][i]>a[pos][i]) pos=j;
      }
      if(a[pos][i]==0){//无解与无数个解都是不可行，求不了逆矩阵
          puts("No Solution");
          return ;
      }
      swap(a[pos],a[i]);
        //同普通高斯消元一样
      int inv=qpow(a[i][i],mod-2);
        //由于要取模，故不可再用小数，于是使用逆元
      for(int j=1;j<=n+n;++j){//注意，不是n+1，而是n+n
          a[i][j]=a[i][j]*inv%mod;
      }
      for(int j=1;j<=n;++j){
          if(j==i) continue;
          int rate=a[j][i];
          for(int k=1;k<=n+n;++k){//同样是n+n，不是n+1
              a[j][k]=(a[j][k]-rate*a[i][k]%mod+mod)%mod;
          }
      }
        //同高斯消元一样
    }
    for(int i=1;i<=n;++i){//输出
        for(int j=1;j<=n;++j){
            write(a[i][j+n]);SPACE;
        }
        endl;
    }
    }

异或方程组

• 介绍

  异或方程组形如$\begin{cases} & a_{1,1}x1\oplus a{1,2}x2\oplus a{1,3}x3\oplus ... \oplus a{1,n}x_n=b1 \ & a{2,1}x1\oplus a{2,2}x2\oplus a{2,3}x3\oplus ... \oplus a{2,n}x_n=b2\ & a{3,1}x1\oplus a{3,2}x2\oplus a{3,3}x3\oplus ... \oplus a{3,n}x_n=b3\ & ...\ & a{m,1}x1\oplus a{m,2}x2\oplus a{m,3}x3\oplus ... \oplus a{m,n}x_n=b_m\ \end{cases}

• 操作

  异或其实就是不进位的加法，我们仍然可以写出增广矩阵，然后使用高斯消元法。

  不同的是，每个变量的取值只有 0 和 1，因此若最终有 q个自由元，方程解的数量为 2^q。

  在编写程序时可以使用二进制来进行状态压缩，或是直接使用 bitset加速

• 代码实现

  注：当然，用高斯消元也是可行的，但是本人偏向喜欢高斯-约旦消元

  void Gauss() {
      int i,j;
      int flag=0;
      //cnt表示枚举的每一列
      //i表示行
      for (i=1;i<=n; i++) {
          int pos=i;
          for (j=i+1;j<=n;j++) { 
              if(a[j][i]){
                  pos=j;
                  break;
              }
              //因为异或方程组的系数无非0或1，所以有1那么一定是最大的
          }
          for(j=1;j<i;j++){
              if(vis[j]&&a[j][i]){
                  pos=j;
                  break;
                  //还需要看前面有些 j 系数都为0 但是 a[j][i] 不一定为0
                  //不能漏掉
              }
          }
          if (!a[pos][i]){//如果该未知数的系数为0，要么无解，要么无数个解(2^p个解) 
              vis[i]=1;
              flag=1;
              continue;//因为剩下的全为0，对于 1~cnt-1 处理不了，于是不处理      
          }
          swap(a[pos],a[i]);//换到最上面
          //系数为1，未必要化系数为1
          for (j=1;j<=n;j++){ //把下面的都变为0
              if(j==i) continue;
              if(a[j][i]){ 
                  for(int k=1;k<=n+1;k++){
                      a[j][k]^=a[i][k];
                  }
              }
          }
      }
      if (flag){ //说明有全为0的情况，即 要么有解，要么有无数个解
          for(i=1;i<=n;i++){
              int f=0;
              for(j=1;j<=n;j++){
                  if(a[i][j]){
                      f=1;
                      break;
                  }
              }
              //先要确定该行 1~n 系数都为0，才能判断 0=？ 的情况
              if(!f&&a[i][n+1]!=0){//先判断无解
                  puts("无解");
                  return ;
              }
          }
          puts("无数个解");//再判断有解
          return ;
      }
      for(i=1;i<=n;i++){
          write(a[i][n+1]);
          //输出解
          endl;
      }
  }

  习题

  线性方程组 逆矩阵 异或方程组 高斯消元法应用