高斯消元

2013-05-12

姿势

数论

数学

线性代数

1.2k 字

3 分钟

高斯消元提供了一种简单的方法解决方程组联立的问题。

基础知识

高斯消元是线性代数中的一类方法，也是一种算法，用来求解线性方程组，并可以求出矩阵的秩，以及求解出可逆方阵的逆矩阵。

根据初等变换原理，将增广矩阵转换为行阶梯矩阵，然后回代求出方程的解。

在程序设计中，常用的算法是： (我们设方程组中方程的个数为equ，变元的个数为var)

把方程组转换成增广矩阵。
利用初等行变换来把增广矩阵转换成行阶梯阵。枚举k从0到equ – 1，当前处理的列为col(初始为0) ，每次找第k行以下(包括第k行)，col列中元素绝对值最大的列与第k行交换。如果col列中的元素全为0，那么则处理col + 1列，k不变。
转换为行阶梯阵，判断解的情况。

① 无解当方程中出现(0, 0, …, 0, a)的形式，且a != 0时，说明是无解的。 ② 唯一解 条件是k = equ，即行阶梯阵形成了严格的上三角阵。利用回代逐一求出解集。 ③ 无穷解 条件是k < equ，即不能形成严格的上三角形，自由变元的个数即为equ – k，但有些题目要求判断哪些变元是不缺定的。

这里单独介绍下这种解法：

首先，自由变元有var - k个，即不确定的变元至少有var - k个。我们先把所有的变元视为不确定的。在每个方程中判断不确定变元的个数，如果大于1个，则该方程无法求解。如果只有1个变元，那么该变元即可求出，即为确定变元。

以上介绍的是求解整数线性方程组的求法，复杂度是O(n^3)。浮点数线性方程组的求法类似，但是要在判断是否为0时，加入EPS，以消除精度问题。

贴代码

代码（雍哥修改的不知道谁的）——整型参数型

#define maxn 405

int equ,var;
//有equ个方程,val个变元 equ行(0~equ-1) val+1列(增广矩阵)(0~val)
int mat[maxn][maxn];//行列式矩阵
int x[maxn];//解
bool free_x[maxn];//是否是不确定变元

void Debug()
{
    for(int i=0;i<equ;++i)
    {
        for(int j=0;j<=var;++j)
            printf("%01d ",mat[i][j]);
        puts("");
    }puts("");
}
inline int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
inline int lcm(int a,int b){return a*b/gcd(a,b);}

//-2 有浮点数解 -1 无解 0有一个解 >0无穷解,返回自由变元个数
int gauss()
{
    int k,max_r,col;//当前行,当前列绝对值最大的行;当前列
    //初始化
    memset(x,0,sizeof(x));
    memset(free_x,true,sizeof(free_x));
    //转换为阶梯矩阵
    for(k=0,col=0;k<equ&&col<var;++k,++col)
    {
        max_r=k;
        for(int i=k+1;i<equ;++i)
            if(abs(mat[i][col])>abs(mat[max_r][col]))max_r=i;
        if(max_r!=k)
            for(int i=k;i<=var;++i)swap(mat[k][i],mat[max_r][i]);
        if(!mat[k][col]){--k;continue;}
        for(int i=k+1;i<equ;++i)if(mat[i][col])
        {
            int l=lcm(abs(mat[i][col]),abs(mat[k][col]));
            int ta=l/abs(mat[i][col]),tb=l/abs(mat[k][col]);
            if(mat[i][col]*mat[k][col]<0)tb=-tb;
            for(int j=col;j<=var;++j)
                mat[i][j]=mat[i][j]*ta-mat[k][j]*tb;
        }
    }//Debug();
    //1.无解的情况 存在行(0,0,...,a) a!=0
    for(int i=k;i<equ;++i)
        if(mat[i][col])return -1;
    //2.无穷解的情况 秩不是equ
    if(k<var)
    {
        for(int i=k-1;i>=0;--i)
        {
            int free_x_num=0,free_index;
            for(int j=0;j<var;++j)
                if(mat[i][j]&&free_x[j])++free_x_num,free_index=j;
            if(free_x_num>1)continue;
            int tmp=mat[i][var];
            for(int j=0;j<var;++j)
                if(mat[i][j]&&j!=free_index)tmp-=mat[i][j]*x[j];
            x[free_index]=tmp/mat[i][free_index];
            free_x[free_index]=false;
        }
        return var-k;
    }
    //3.唯一解 严格三角阵 满秩
    for(int i=var-1;i>=0;--i)
    {
        int tmp=mat[i][var];
        for(int j=i+1;j<var;++j)
            if(mat[i][j])tmp-=mat[i][j]*x[j];
        if(tmp%mat[i][i])return -2;//有浮点数
        x[i]=tmp/mat[i][i];
    }
    return 0;
}

对于浮点参数型，基本实现形式差不多，只不过求解的时候注意点浮点数的运算就行了。

代码（来自这里）：

#define MAXN 405
const double EPS=1e-9;

double A[MAXN][MAXN];
double X[MAXN];
bool Free[MAXN];
int equation, variable;

void debug()
{
    printf("The Matrix is :\n");
    for(int i=0; i<=variable+1; i++)
    {
        printf("%8d",i);
        if(i==0) printf("|");
    }
    printf("\n");
    for(int i=0; i<=variable+1; i++)
    {
        printf("-------");
        if(i==0) printf("-+");
        else printf("-");
    }
    printf("\n");

    for(int i=1; i<=equation; i++)
    {
        printf("%8d|",i);
        for(int j=1; j<=variable+1; j++)
        {
            printf("%8.2lf",A[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

/***************************************************************************/
///   功能说明：
///   1、不需带入参数，参数已经在 build() 函数中设置好了
///   2、函数返回 -1 ：无解
///   3、函数返回 0  ：有且仅有一组解
///   4、函数返回 其他正数 ：有多解，返回的值是不确定变元的个数
/**************************************************************************/
int gauss()
{
    int col=1, row=1;
    for(; row<=equation && col<=variable; row++, col++)
    {
        int max_row=row;

        //找到 col 那列元素绝对值最大的那行与当前行交换，减小精度误差
        for(int i=row+1; i<=equation; i++)
            if(abs(A[i][col])>abs(A[max_row][col]))
                max_row=i;

        //如果 col 那列元素最大是 0，表明这一列全部是 0，处理下一列
        if(abs(A[max_row][col])<=EPS)
        {
            row--;
            continue;
        }

        //如果不是同一行，交换元素
        if(max_row!=row)
        {
            for(int i=col; i<=variable+1; i++)
                swap(A[max_row][i], A[row][i]);
        }

        //枚举要删去的行
        //for(int i=row+1; i<=equation; i++)
        for(int i=1; i<=equation; i++)
        {
            if(i==row) continue;
            if(abs(A[i][col])<=EPS) continue;

            double ta=A[row][col];
            double tb=A[i][col];

            for(int j=1; j<=variable+1; j++)
                A[i][j]=A[i][j]*ta-A[row][j]*tb;
        }
//        printf("row=%d col=%d\n", row, col);
//        debug();
    }

//    debug();

    // 1、没有解的情况：
    for(int i=row; i<=equation; i++)
        if(A[i][variable+1]!=0) return -1; //表明无解

    // 2、无穷解的情况：
    if(row<variable+1)
    {
        for(int i=row-1; i>=1; i--)
        {
            int free_x_num=0, free_index=-1;
            for(int j=1; j<=variable; j++)
                if(Free[j] && A[i][j]!=0)
                {
                    free_x_num++;
                    free_index=j;
                }
            if(free_x_num>1) continue;

            Free[free_index]=false;

            double temp=A[i][variable+1];
            for(int j=1; j<=variable; j++)
                if(!(abs(A[i][j])<=EPS) && j!=free_index)
                    temp-=A[i][j]*X[j];
            X[free_index]=temp/A[i][free_index];
        }
        return variable-row+1;
    }

    // 3、唯一解的情况：
    for(int i=variable; i>=1; i--)
    {
        double temp=A[i][variable+1];
        for(int j=i+1; j<=variable; j++)
            if(!(abs(A[i][j])<=EPS))
                temp-=A[i][j]*X[j];
        X[i]=temp/A[i][i];
    }
    return 0;
}

void build(int n, int m)
{
    memset(Free, true, sizeof(Free));
    memset(X, 0, sizeof(X));
    variable=n, equation=m;
    for(int i=1; i<=equation; i++)
        for(int j=1; j<=variable+1; j++)
            scanf("%lf",&A[i][j]);
    debug();
}

void PRIMAXNT(int flag)
{
    if(flag==-1)
    {
        printf("MAXNO Solution !\n");
        return;
    }
    if(flag==0)
    {
        for(int i=1; i<=variable; i++)
            printf("x%d : %4.4lf\n", i, X[i]);
        return;
    }
    printf("There are %d variables who are indeterminate !\n",flag);
    for(int i=1; i<=variable; i++)
    {
        printf("x%d : ", i);
        if(Free[i]) printf("indeterminate !\n");
        else printf("%4.4lf\n",X[i]);
    }
}

主要解决的问题

高斯消元主要解决这几类问题：

普通的消元问题

给出几个方程组，求出每个解。一般套用模板就可以了。

染色问题

一个方格只有01两种状态，改变其中一个方格会影响其它周围方格的状态。

对于这种问题，首先要建立方程组。关于方程组的建立，首先要把方格转换成一行，然后构造一个 n * n 的矩阵。那么对于矩阵中的( i,j )就表示j可以影响i的状态。

那么如果 i 能影响 j 点，就应该把( j,i )设为1。

用高斯消元的时候，注意把中间过程 mod 2，因为只有 01 两种状态。

mod p问题

方程组对p取余，相当于一个线性同余型的问题吧。

对于这种问题，分为两种情况：

如果p为质数，这样就可以参考染色问题，根据性质，直接将中间过程 mod p，最后得出的解保证小于 p 即可。
如果p为 非质数 ，这样就不能对中间过程 mod p，因为这样会形成循环节而扩大解集。我也没有什么好的解决方法，暂时可以中间过程不要进行任何取余操作。求解集的时候可以枚举 0~p 的数字，满足同余关系就可以作为解的一部分。这个也要注意解不要超过 p 。

找出逆矩阵

高斯消元法可以用来找出一个可逆矩阵的逆矩阵。设 A 为一个 N * N 的矩阵，其逆矩阵可被两个分块矩阵表示出来。将一个 N * N 单位矩阵放在 A 的右手边，形成一个 N * 2N 的分块矩阵 B = [ A,I ] 。经过高斯消元法的计算程序后，矩阵 B 的左手边会变成一个单位矩阵 I ，而逆矩阵 A - 1 会出现在B 的右手边。

假如高斯消元法不能将A 化为三角形的格式，那就代表A 是一个不可逆的矩阵。

应用上，高斯消元法极少被用来求出逆矩阵。高斯消元法通常只为线性方程组求解。

#高斯消元

#基础知识

#贴代码

#主要解决的问题

#普通的消元问题

#染色问题

#mod p问题

#找出逆矩阵

#参考内容