主题

主题

 
   

从高斯消元说起

一些方便处理的概念

向量:我们可以把它理解为一个数组。

比方说在下图中,x是一个向量,y是一个向量。至于a的话,我们用矩阵来称呼好了。

方阵:我们可以把它理解为一个长等于宽的矩阵。

关于矩阵的下标:可能有的时候用0开始有点时候用1开始不必介意。只需要大致明白此意即可。在实际运用中,这是细节问题。

乘法:我们把向量视为矩阵,可以把向量之间的乘法视为矩阵乘法。

主对角线:我们可以视为是左上-右下对角线。


  • 高斯消元

一个处理线性方程组的方法。

用增广矩阵和两个向量来实现,类似于:


这是一般的表现形式。

给出左边的M方阵,以及右边的y向量,求中间的x向量。

一般的做法就是消元成一个上三角矩阵。

搬一个代码:

void Gauss(func &M,double *y,double *x,int len) {

    int vi,vj,vk;

    for(vk=0; vk<len; vk++) {

        for(vi=vk+1; vi<len; vi++) {

            M.A[vi][vk]/=M.A[vk][vk];

            for(vj=vk+1; vj<len; vj++)M.A[vi][vj]-=M.A[vi][vk]*M.A[vk][vj];

            y[vi]-=M.A[vi][vk]*y[vk];

        }

    }

    for(vi=len-1; vi>=0; vi--) {

        x[vi]=y[vi];

        for(vj=vi+1; vj<len; vj++)x[vi]-=x[vj]*M.A[vi][vj];

        x[vi]/=M.A[vi][vi];

    }

}

这个做法比较大众化。

有的时候我们会采用列主元的方法。这既可以避免较大的精度误差,也可以处理掉首项为0的情况。

也就是把首项绝对值最大行与当前行交换。

在高斯消元中,对解不会造成任何影响。但毕竟交换了矩阵的行,在后面的讨论中,我们会涉及到它的影响。

再搬一个代码:

void Gauss(func &M,double *x,double *y,const int &len) {

    int vi,vj,vk,pos;

    double maxn;

    for(vk=1; vk<=len; vk++) {

        for(maxn=0,vi=vk; vi<=len; vi++)if(maxn<std::abs(M.A[vi][vk])) {

                maxn=std::abs(M.A[vi][vk]);

                pos=vi;

            }

        for(vi=vk; vi<=len; vi++)std::swap(M.A[vk][vi],M.A[pos][vi]);//这里只交换了有用信息。在之后的讨论中,可能需要将整行交换

        std::swap(y[vk],y[pos]);//注意要交换y

        for(vi=vk+1; vi<=len; vi++) {

            M.A[vi][vk]/=M.A[vk][vk];

            for(vj=vk+1; vj<=len; vj++)M.A[vi][vj]-=M.A[vi][vk]*M.A[vk][vj];

            y[vi]-=M.A[vi][vk]*y[vk];

        }

    }

    for(vi=len; vi>=1; vi--) {

        x[vi]=y[vi];

        for(vj=vi+1; vj<=len; vj++)x[vi]-=x[vj]*M.A[vi][vj];

        x[vi]/=M.A[vi][vi];

    }

}

在接下来的讨论中,一般M矩阵被视为一个方阵。


  • 上三角矩阵与下三角矩阵

很容易理解。

上三角矩阵指的是:一个方阵中对角线以及对角线以上有元素的特殊的“矩阵”。

下三角矩阵指的是:一个方阵中对角线以及对角线以下有元素的特殊的“矩阵”。

注意,这里的对角线是主对角线。

解线性方程组所使用的高斯消元本质上只是在方程组中大力消元而已,并不能挖掘出矩阵的各种性质。

接下来我们从方阵上做一些研究。


  • LU分解

我们设A=LU,称为A的LU分解。

其中L是一个单位下三角矩阵(即主对角线为1,主对角线下可能有非0的元素),U是一个上三角矩阵。

y=Ax=LUx。x与y是向量。比方说:


(注意U矩阵的右下角是可以有元素只是值为0)

设z=Ux。那么我们可以得到两个式子

  1. z=Ux

  2. y=Lz

因为U是一个上三角矩阵,所以如果能求出z那么我们可以很快地求出x。


这就是高斯消元的第二部分。

因为L是一个单位下三角矩阵,所以我们可以更快地求出来z。


用眼睛都可以看出来。

比较重要的发现是,这两个部分都是O(n^2)的。

也就是意味着,我们要是求出了LU分解,就可以在O(n^2)的时间内,对于一个向量y,求出对应的x。

那么怎么求解这个LU分解呢?

幸运的是,前人已经给出了求解矩阵LU分解的方法,而且这个方法很简单。

更幸运的是,求解LU分解的方法和之前的高斯消元是一样的。

在做完求解线性方程组的高斯消元之后,A数组的对角线及对角线上就是U矩阵,对角线下就是L矩阵。

搬一个代码:

void LU(func &M,const int &len) {

    int vi,vj,vk;

    for(vk=1; vk<=len; vk++)

        for(vi=vk+1; vi<=len; vi++) {

            M.A[vi][vk]/=M.A[vk][vk];

            for(vj=vk+1; vj<=len; vj++)M.A[vi][vj]-=M.A[vi][vk]*M.A[vk][vj];

        }

}

在没有给出x与y的时候做LU分解。代码唯一的不同就是后面没有处理y向量而已。

拿到了LU分解之后,自然是要解方程。

先解出z向量,再解出x向量。

搬一个代码:

double z[210];

void Gauss(const func &M,double *x,double *y,const int &len) {

    int vi,vj;

    for(vi=1; vi<=len; vi++) {

        z[vi]=y[vi];

        for(vj=1; vj<=vi-1; vj++)z[vi]-=z[vj]*M.A[vi][vj];

    }

    for(vi=len; vi>=1; vi--) {

        x[vi]=z[vi];

        for(vj=vi+1; vj<=len; vj++)x[vi]-=x[vj]*M.A[vi][vj];

        x[vi]/=M.A[vi][vi];

    }

}

已经很明显地看出来了,两部分代码合并基本上就是原来的高斯消元。LU分解只是借助了一个z向量,把之前O(n^3)的部分提出来,变为了预处理。这样解多次高斯消元就变成了了O(n^3+kn^2)的了。


  • LUP分解

注意要记录P,来处理y的变化。

也就是说,我们仍然采取列主元的方法。这样可能会交换行。我们用一个P数组(其实是一个置换矩阵)记录一下这个当前行原来是第几行。

如果需要那LUP分解来做多次的高斯消元的话,就需要记录一下P,在初次进入消元的时候把y按照P交换一下。

搬一个代码:

void LUP(func &A,const int &len) {

    int vi,vj,vk,pos;

    double maxn;

    for(vi=1; vi<=len; vi++)P[vi]=vi;

    for(vk=1; vk<=len; vk++) {

        for(maxn=0,vi=vk; vi<=len; vi++)if(maxn<std::abs(A.A[vi][vk])) {

                pos=vi;

                maxn=std::abs(A.A[vi][vk]);

            }

        std::swap(P[vk],P[pos]);

        for(vi=1; vi<=len; vi++)std::swap(A.A[vk][vi],A.A[pos][vi]);//注意这里要全部交换。切记切记。

        for(vi=vk+1; vi<=len; vi++) {

            A.A[vi][vk]/=A.A[vk][vk];

            for(vj=vk+1; vj<=len; vj++)A.A[vi][vj]-=A.A[vi][vk]*A.A[vk][vj];

        }

    }

}

void Gauss(func &A,double *x,double *y,const int &len) {

    int vi,vj;

    for(vi=1; vi<=len; vi++)tmp[vi]=y[P[vi]];//tmp数组记录变化后的y数组

    for(vi=1; vi<=len; vi++)y[vi]=tmp[vi];

    for(vi=1; vi<=len; vi++) {

        z[vi]=y[vi];

        for(vj=1; vj<=vi-1; vj++)z[vi]-=z[vj]*A.A[vi][vj];

    }

    for(vi=len; vi>=1; vi--) {

        x[vi]=z[vi];

        for(vj=vi+1; vj<=len; vj++)x[vi]-=x[vj]*A.A[vi][vj];

        x[vi]/=A.A[vi][vi];

    }

}

交换线性方程组是不会影响x的顺序的,所以不用交换x。


(这一部分的模板题可以在luogu提交:https://www.luogu.org/problem/show?pid=3389#sub


  • 行列式

一个O(n!)的行列式计算方法:

对于一个方阵,按照字典序列出长度为n的所有的排列。对于每一个排列,若A[i]=j,则取方阵第i行第j列的元素。记这个排列的sum为所有取出的元素的积。

如果这个排列是第奇数项排列,则令行列式+=sum,反则令行列式-=sum。

比方说对于一个n=3的矩阵:



  • LU分解计算行列式

用到两个性质(定理):

  1. 对于方阵A,X,Y,如果A=XY,那么del(A)=del(X)*del(Y)(del即为行列式的值)。(并不会证明,也许以后会补上)

  2. 对于一个三角矩阵A,del(A)为对角线上所有元素的乘积。

那么我们得到了LU分解之后,就可以直接计算出当前方阵的行列式了。

或者更加直接的说,若A=LU,则del(A)=del(U)=pai(U[i][i])。

如何处理模意义下的行列式?

辗转相除法。

搬一下代码:

long long sum=1;

void LU(func &M,const int &len) {

    int vi,vj,vk,tmp;

    for(vi=1; vi<=len; vi++)

        for(vj=1; vj<=len; vj++)(M.A[vi][vj]+=intp)%=intp;

    for(vk=1; vk<=len; vk++) {

        if(!M.A[vk][vk]) {

            for(vi=vk+1; vi<=len; vi++)if(M.A[vi][vk])break;

            for(vj=vk; vj<=len; vj++)std::swap(M.A[vi][vj],M.A[vk][vj]);

            (sum*=(intp-1))%=intp;

        }

        for(vi=vk+1; vi<=len; vi++) {

            while(M.A[vi][vk]) {

                tmp=M.A[vi][vk]/M.A[vk][vk];

                for(vj=vk; vj<=len; vj++)(((M.A[vi][vj]-=tmp*M.A[vk][vj])%=intp)+=intp)%=intp;

                if(!M.A[vi][vk])break;//注意如果已经达到目的了就不交换了,行列式也不用再乘-1了。

                (sum*=(intp-1))%=intp;

                for(vj=vk; vj<=len; vj++)std::swap(M.A[vi][vj],M.A[vk][vj]);

            }

        }

        (sum*=M.A[vk][vk])%=intp;

    }

}

为什么这么长呢?因为模意义下的行列式求值,不像高斯消元解浮点数的线性方程组,这个问题在实践中基本上会出现都当前行首项为0的情况,需要把行交换。

好在行列式满足了两个性质,让我们可以使用辗转相除以及行交换:

  1. 互换行列的任意两行(两列)行列式变号。

  2. 把行列式的任一行(列)的元素乘以同一个数后,加到另一行(列)的互对应元素上去,行列式不变。

注意,变号——即乘以-1,事实上应该是乘以(intp-1)。

当然,在这里使用列主元的方法也是可以的。注意如果当前行的首项就是主元,不要乘以-1。

事实上它是LUP分解……只是这里P不需要记录。


(这一部分的模板题可以在SPOJ上提交:http://www.spoj.com/problems/DETER3/


  • 矩阵-树定理

如何求图的生成树的个数?

结论:

令E为G的邻接矩阵(元素代表边的个数,可以有重边与自环),D为度数矩阵(对角线为度数),M=D-E。

M去掉任何一行一列之后的行列式的值就是图的生成树的个数。

然后我们就可以O(n^3)求出图的生成树的个数了。

刚开始以为需要绝对值,直到后来看懂了巨佬的博客,其实是可以不用取绝对值的。因此一些模意义下的行列式问题就变得方便了。

分享:http://blog.csdn.net/u013010295/article/details/47451451(尊重作者,不引用以及转载了)


(这一部分的模板题可以在SPOJ上提交:http://www.spoj.com/problems/HIGH/


  • 矩阵求逆

一般用到的是对方阵的求逆。这个也是O(n^3)的。在我们学习了LU分解之后。

然而对矩阵的求逆也是基于高斯消元。

如果是用于矩阵收敛这一部分的话,其实作用也不大,因为一般这样的题目会做一个矩阵的变化变为高斯消元。不需要用到矩阵求逆。

但是矩阵求逆的范围自然更广。

先疏通一下定义:

对于一个A矩阵,如果存在A*B=1,那么B为A的逆矩阵。

虽然矩阵乘法不满足交换律,但是我们发现A与B是互为逆矩阵的。其实这是可以用结合律来解释:

已知A*B=1,

那么:

A*(B*A)=(A*B)*A=1*A=A

所以:

(B*A)*(B*A)=B*(A*(B*A))=(B*A)

等式连边同时乘以(B*A)的逆矩阵,即可得到:

B*A=1

一个朴素的想法是把逆矩阵的所有元都列出方程,这样做的复杂度是O(n^6)的。

事实上有更加高效的方法,请看下例:


假设我们现在求的是A方阵的逆方阵B。其实拆开来看,就可以发现:


本质上是在高斯消元。

由于我们要对同一个方阵做多次高斯消元,所以我们要用到LU分解。

所以时间复杂度仍然为O(n^3)。

搬一个代码:

void LUP(func &A,const int &len) {

    int vi,vj,vk,pos;

    double maxn;

    for(vi=1; vi<=len; vi++)P[vi]=vi;

    for(vk=1; vk<=len; vk++) {

        for(maxn=0,vi=vk; vi<=len; vi++)if(maxn<std::abs(A.A[vi][vk])) {

                pos=vi;

                maxn=std::abs(A.A[vi][vk]);

            }

        std::swap(P[vk],P[pos]);

        for(vi=1; vi<=len; vi++)std::swap(A.A[vk][vi],A.A[pos][vi]);

        for(vi=vk+1; vi<=len; vi++) {

            A.A[vi][vk]/=A.A[vk][vk];

            for(vj=vk+1; vj<=len; vj++)A.A[vi][vj]-=A.A[vi][vk]*A.A[vk][vj];

        }

    }

}

double x[120],y[120],z[120],tmp[120];

func tmpA;

void Gauss(func &A,double *x,double *y,const int &len) {

    int vi,vj;

    for(vi=1; vi<=len; vi++)tmp[vi]=y[P[vi]];

    for(vi=1; vi<=len; vi++)y[vi]=tmp[vi];

    for(vi=1; vi<=len; vi++) {

        z[vi]=y[vi];

        for(vj=1; vj<=vi-1; vj++)z[vi]-=z[vj]*A.A[vi][vj];

    }

    for(vi=len; vi>=1; vi--) {

        x[vi]=z[vi];

        for(vj=vi+1; vj<=len; vj++)x[vi]-=x[vj]*A.A[vi][vj];

        x[vi]/=A.A[vi][vi];

    }

}

void inver(const func &A,func &res,const int &len) {

    int vi,vj;

    tmpA=A;//这里用for自然会好一点

    LUP(tmpA,len);

    for(vi=1; vi<=len; vi++) {

        for(vj=1; vj<=len; vj++)y[vj]=0;

        y[vi]=1;

        Gauss(tmpA,x,y,len);

        for(vj=1; vj<=len; vj++)res.A[vj][vi]=x[vj];

    }

}

为了不改变A,所以就新令了一个tmpA。


我们得到经验:

  1. 高斯消元可以不列主元但是一般列主元。

  2. 求行列式必须要列主元(不列主元也要判断首项是否为0,主要是处理模意义下的行列式的求值)。

  3. 矩阵求逆可以不列主元(LU分解)但是一般列主元(LUP分解)。

  4. 在高斯消元/计算行列式的时候可以只从vk开始交换。但是,在LUP分解的时候一定要把整个行都置换。交换对三角矩阵有很大影响。在不能调试纯靠熟练的的考场,这很有可能就是100到0分的罪魁祸首!切记切记。

 
 
评论
 
 
热度(2)