高斯消元学习笔记

概述

高斯消元用于求解线性方程组，也可用于求矩阵的秩和逆。
方法是先将增广矩阵消为上三角行列式，再逐一回代求解。
高斯-约旦消元法：在消元过程中直接消为最简行阶梯型，省略了回代的过程。
下文中高斯消元的算法与代码实际上都是高斯-约旦消元法。

注意：
1. 当系数矩阵的秩小于增广矩阵的秩时，方程组无解
2. 当系数矩阵的秩等于增广矩阵的秩等于变量数时，方程组有唯一解。
3. 当系数矩阵的秩等于增广矩阵的秩小于变量数时，方程组有无穷解。

算法

高斯消元函数接受一个增广矩阵，然后返回解的存在情况，以及有唯一解时，获得解的方式。
$vars$：变量数
$equs$：方程数
$frees$：自由变元数
1. 获取初始$frees$。若$vars>equs$，不可能有唯一解，此时$frees=vars-equs$，否则为$0$.
2. 逐行消元，需要进行$min(vars,equs)$次，记录当前行为$row$.
2.1 找到目标行。从$row$开始遍历所有的行，找到第$row$列最大的一个，并将该行与当前行交换。
2.2 如果目标行的第$row$列为$0$（即所有剩余行第$row$列都为$0$）且常数列不为$0$，返回无解。
2.3 如果目标行的第$row$列为$0$且常数列也为$0$，$frees+=1$，开始消下一行。
2.4 对于除目标行之外的所有行，每列均减去第$row$列与目标行第$row$列的比值乘上目标行当前列的比值，即将第$row$列消为$0$.
3. 若$equs>vars$，仍需要判断剩余的等式，若常数项不为$0$，返回无解。
4. 返回$frees$，此时若$frees$为$0$表示有唯一解，且每行的常数项除以对应列即为该项的解。若$frees$不为$0$，即表示自由变元的个数。

不同种类的高斯消元代码及例题

浮点高斯消元

最常规的高斯消元

//Gauss_float.cpp 浮点高斯消元
/**
* 浮点高斯消元,传入增广矩阵,求解线性方程组
* 返回值:-1表示无解,0表示唯一解,大于0表示有这么多自由元
* 若有唯一解,mat[i][vars]/mat[i][i]表示第i个解
* 步骤:交换到当前列最大的行,判断是否为0,乘除消元
*/
int Gauss(vector<vector<double>> &mat)
{
    const int equs = mat.size(), vars = mat[0].size() - 1;
    int frees = max(0, vars - equs);

    for(int row = 0; row < min(vars,equs); row++)
    {
        int mx = row;
        for(int trow = row + 1; trow < equs; trow++)
            if(fabs(mat[trow][row]) > fabs(mat[mx][row]))
                mx = trow;
        if(mx != row) swap(mat[row], mat[mx]);

        if(fabs(mat[row][row]) < eps) 
        {
            if(fabs(mat[row][vars]) > eps) return -1;
            frees++;
            continue;
        }

        for(int trow = 0; trow < equs; trow++)
        {
            if(trow == row) continue;
            double mul = mat[trow][row] / mat[row][row];
            for(int col = row; col <= vars; col++)
                mat[trow][col] -= mat[row][col] * mul;
        }
    }
    for(int row = vars; row < equs; row++)
        if(fabs(mat[row][vars]) > eps)
            return -1;

    return frees;
}

luogu P3389

模板题

二进制高斯消元

操作全部改成异或，找目标行时只要找到当前列为1的即可

/**
* 二进制高斯消元,传入增广矩阵,求解异或方程组
* 返回值:-1表示无解,0表示唯一解,大于0表示有这么多自由元
* 若有唯一解,mat[i][vars]表示第i个解
*/
int Gauss(vector<vector<bool>> &mat)
{
    const int equs = mat.size(), vars = mat[0].size() - 1;
    int frees = max(0, vars - equs);

    for(int row = 0; row < min(vars,equs); row++)
    {
        int flag = 0;
        for(int trow = row; trow < equs; trow++) if(mat[trow][row])
        {
            if(trow != row) swap(mat[row], mat[trow]);
            flag = 1;
            break;
        }
        if(flag==0) 
        {
            if(mat[row][vars]) return -1;
            frees++;
            continue;
        }
        for(int trow = 0; trow < equs; trow++)
        {
            if(trow == row || !mat[trow][row]) continue;
            for(int col = row; col <= vars; col++)
                mat[trow][col] = mat[row][col] ^ mat[trow][col];
        }
    }
    for(int row = vars; row < equs; row++)
        if(mat[row][vars])
            return -1;

    return frees;
}

POJ 1222

5*6=30个灯，按一个开关可以将自身及上下左右最多四个灯转换亮暗，求如何能将给定的亮暗状态转变为全暗。

列方程调用函数即可，实际上每个开关只会按1次或0次，而且上一行可以完全确定下一行，所以可以通过枚举第一行来确定整个状态。

const int M = 30;
int main(void)
{
    vector<vector<bool>> save(M, vector<bool>(M + 1));
    const int go[5][2] = {{0,0},{0,1},{0,-1},{1,0},{-1,0}};
    for(int pos = 0; pos < M; pos++)
    {
        int x = pos / 6, y = pos % 6;
        for(int k = 0; k < 5; k++)
        {
            int nx = x + go[k][0], ny = y + go[k][1];
            if(nx >= 0 && nx < 5 && ny >= 0 && ny < 6)
                save[pos][nx * 6 + ny] = 1;
        }
    }

    int T=read();
    for(int t=1;t<=T;t++)
    {
        vector<vector<bool>> mat = save;
        for(int i = 0; i < M; i++)
            mat[i][M] = read();
        Gauss(mat);
        printf("PUZZLE #%d\n",t );
        for(int i=0;i<M;i++)
            printf("%d%c",(int)mat[i][M],i%6==5?'\n':' ' );

    }

    return 0;
}

POJ 1830

求可行的方法数，实际上就是自由变元个数的任意组合。

int main(void)
{
    int T = read(), ta, tb;
    while(T--)
    {
        int n = read();
        vector<vector<bool>> mat(n, vector<bool>(n + 1));
        for(int i = 0; i < n; i++)
        {
            mat[i][n] = read();
            mat[i][i] = 1;
        }
        for(int i = 0; i < n; i++)
            mat[i][n] = mat[i][n] ^ read();
        while(scanf("%d%d", &ta, &tb) && ta)
            mat[tb - 1][ta - 1] = 1;
        int res = Gauss(mat);
        if(res == -1)
            printf("Oh,it's impossible~!!\n");
        else
            printf("%d\n", 1 << res );
    }

    return 0;
}

综合应用

HDU 4418 Time travel 期望dp,高斯消元

n(100)个点的线段，标号为0到n-1,每次有pi的几率走i步(1<=i<=m(100)),走到头时转向,给定起点和终点,以及起点时的朝向,问停在终点的步数期望?

需要注意的点:
1. 转向不好处理,n个点有n-2种状态,不如扩展成n-2个点,走到尾后从头开始走.
扩展公式:y=2n-2-x,因为起点给定了方向,仍为一个点.但终点可以有两个,求期望时应当计算平均值.
接下来是期望dp:
状态表示:dp[x]表示从起点走到xi停下的步数期望.
状态转移:$dp_x = \Sigma_{i=1}^{m}p_i*(dp_{x+i} + i)$,注意全正走和全倒走是一样的.