「算法笔记」高斯消元

引子

高斯消元$(gauss\ elimination)$是线性代数中的一个算法。可用来求解线性方程组，并可以求出矩阵的秩，以及求出可逆方阵的逆矩阵。
高斯消元法的原理是：用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解。

---

例题

求解一个$n$元线性方程组。

---

引导

想一想我们是如何求解二元一次方程组的。

$\{^{ax+by=c\ (1)}_{dx+ey=f\ (2)}$

保证$a,b,c,d,e,f>0$且$ae\ne bd$

$\ \ \ \ \ \ \ (1)-\frac{a}{d}\times(2)$

$=>(b-\frac{ae}{d})\ y=c-\frac{af}{d}$

$=>y=\frac{cd-af}{bd-ae}$

代入$(1)$得$x=\frac{ce-bf}{bd-ae}$

消元，然后带入，对不对？

现在我们就要使用编程模拟这个过程。

---

算法流程

方程组为：
$a_{1,1}x_1+a_{1,2}x_2+...+a_{1,n}x_n=b_1$
$a_{2,1}x_1+a_{2,2}x_2+...+a_{2,n}x_n=b_2$
$...$
$a_{n,1}x_1+a_{n,2}x_2+...+a_{n,n}x_n=b_n$

• 进行第$i$次消元时，先取出任意一行，记为$row$。
• 然后，对于所有没有被取出的行$j$，对其进行消元。具体过程如下：
  
  • 令$d=\frac{a_{j,i}}{a_{row,i}}$。
  • $b_j\leftarrow b_j-d\times b_i$。
  • $a_{j,k}\leftarrow a_{j,k}-d\times a_{row,k}$。
  • 等到消元完毕时，将$x_{n,n-1,...,1}$依次代入，即可求出方程的解。
    算法时间复杂度$O(n^3)$

  • ---

    优化

    每一次找一个使得$|a_{row,i}|$最大的$row$，可以有效减小精度误差。

    ---

    代码

    //  高斯消元解 n 元线性方程组
    #include <cmath>
    #include <cstdio>
    #include <iostream>
    const double eps = 1e-18;
    const int maxn = 105;
    using namespace std;
    typedef double db;
    typedef long double ldb;
    db t;
    int n;
    ldb a[maxn][maxn], b[maxn], x[maxn];
    int main() {
        scanf("%d", &n);
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                scanf("%lf", &t);
                a[i][j] = t;
            }
            scanf("%lf", &t);
            b[i] = t;
        }
        for (int i = 1; i <= n; i++) {
            int ind = i;
            for (int j = i + 1; j <= n; j++) {
                if (abs(a[j][i]) > abs(a[ind][i])) {
                    ind = i;
                }
            }
            if (abs(a[ind][i]) < eps) {
                puts("-1");
                return 0;
            }
            if (ind != i) {
                swap(b[i], b[ind]);
                for (int j = 1; j <= n; j++) {
                    swap(a[i][j], a[ind][j]);
                }
            }
            for (int j = i + 1; j <= n; j++) {
                double d = a[j][i] / a[i][i];
                b[j] -= b[i] * d;
                for (int k = i; k <= n; k++) {
                    a[j][k] -= a[i][k] * d;
                }
            }
            /*
            for (int j = 1; j <= n; j++) {
                for (int k = 1; k <= n; k++) {
                    printf("%lf ", a[j][k]);
                }
                printf("%lf\n", b[j]);
            }
            */
        }
        for (int i = n; i > 0; i--) {
            for (int j = i + 1; j <= n; j++) {
                b[i] -= x[j] * a[i][j];
            }
            x[i] = b[i] / a[i][i];
        }
        for (int i = 1; i <= n; i++) {
            t = x[i];
            printf("x[%d] = %.6lf;\n", i, t);
        }
        return 0;
    }

    ---

    练习

    求$n\times n$的矩阵的行列式。

    ---

    提示

    只要将矩阵转化成上三角矩阵即可。

    ---

    代码

    //  高斯消元求矩阵行列式
    #include <cmath>
    #include <cstdio>
    #include <iostream>
    const int maxn = 105;
    const double eps = 1e-6;
    typedef long double ldb;
    using namespace std;
    double t;
    int n;
    ldb a[maxn][maxn];
    double read() {
        scanf("%lf", &t);
        return t;
    }
    void write(ldb x) {
        t = x;
        printf("%lf;\n", t);
    }
    int main() {
        scanf("%d", &n);
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                a[i][j] = read();
            }
        }
        ldb x = 1;
        for (int i = 1; i <= n; i++) {
            //  为了减小精度误差，取最大的系数
            int ind = i;
            for (int j = i + 1; j <= n; j++) {
                if (abs(a[ind][i]) < abs(a[j][i])) {
                    ind = j;
                }
            }
            if (ind != i) {
                x = -x;
                for (int j = 1; j <= n; j++) {
                    swap(a[i][j], a[ind][j]);
                }
            }
            if (abs(a[i][i]) < eps) {
                puts("0");
                return 0;
            }
            for (int j = i + 1; j <= n; j++) {
                ldb d = a[j][i] / a[i][i];
                for (int k = i; k <= n; k++) {
                    a[j][k] -= d * a[i][k];
                }
            }
            printf("Step #%d of Gaussian Elimination\n", i);
            for (int j = 1; j <= n; j++) {
                for (int k = 1; k <= n; k++) {
                    t = a[j][k];
                    printf("%lf%c", t, k == n ? '\n' : ' ');
                }
            }
        }
        for (int i = 1; i <= n; i++) {
            x *= a[i][i];
        }
        printf("Det(A) = ");
        write(x);
        return 0;
    }

    ---

    高斯消元与动态规划

    例题

    【LightOJ1511】Snakes and Ladders

    题解

    如果没有蛇，概率动态规划就无后效性，推一下状态转移方程即可。
    如果有蛇，我们就视$dp$值为未知数，列出一个$n$元的线性方程组，用高斯消元法来求解即可。

    代码

    #include <cmath>
    #include <cstdio>
    #include <cstring>
    #include <iostream>
    typedef long double ldb;
    const ldb eps = 1e-6;
    const int maxn = 105;
    const int n = 100;
    using namespace std;
    int tasks, m, to[maxn];
    ldb a[maxn][maxn], b[maxn], x[maxn];
    int main() {
        scanf("%d", &tasks);
        for (int task = 1; task <= tasks; task++) {
            memset(to, 0, sizeof(to));
            for (int i = 1; i <= n; i++) {
                for (int j = 1; j <= n; j++) {
                    a[i][j] = 0;
                }
                b[i] = 0, x[i] = 0;
            }
            scanf("%d", &m);
            for (int u, v, i = 1; i <= m; i++) {
                scanf("%d %d", &u, &v);
                to[u] = v;
            }
            for (int i = 1; i <= n; i++) {
                if (i == n) {
                    a[i][i] = 1;
                } else if (to[i]) {
                    a[i][i] = 1, a[i][to[i]] = -1;
                } else if (i + 6 <= 100) {
                    for (int j = i + 1; j <= i + 6; j++) {
                        a[i][j] = -1;
                    }
                    a[i][i] = 6, b[i] = 6;
                } else {
                    a[i][i] = 100 - i;
                    for (int j = i + 1; j <= n; j++) {
                        a[i][j] = -1;
                    }
                    b[i] = 6;
                }
            }
            for (int i = 1; i <= n; i++) {
                int ind = i;
                for (int j = i + 1; j <= n; j++) {
                    if (abs(a[j][i]) > abs(a[ind][i])) {
                        ind = j;
                    }
                }
                if (i != ind) {
                    swap(b[i], b[ind]);
                    for (int j = 1; j <= n; j++) {
                        swap(a[i][j], a[ind][j]);
                    }
                }
                for (int j = i + 1; j <= n; j++) {
                    ldb d = a[j][i] / a[i][i];
                    b[j] -= b[i] * d;
                    for (int k = 1; k <= n; k++) {
                        a[j][k] -= a[i][k] * d;
                    }
                }
            }
            double tmp;
            for (int i = n; i >= 1; i--) {
                for (int j = i + 1; j <= n; j++) {
                    b[i] -= a[i][j] * x[j];
                }
                x[i] = b[i] / a[i][i];
                tmp = x[i];
            }
            tmp = x[1];
            printf("Case %d: %.10lf\n", task, tmp);
        }
        return 0;
    }
    /*
    100
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    */

    P.S.

    其实，这种类型的题目有时也可以用迭代的方法来做，最后往往也能达到题目要求的精度。

    ---

    总结

    碰到一些有后向性的概率动态规划$(n\le 500)$，首先就应该想到高斯消元。这是信息学竞赛中的常用手段。