zoukankan      html  css  js  c++  java
  • 【Learning】矩阵树定理 Matrix-Tree

    矩阵树定理 Matrix Tree

      
    ​  矩阵树定理主要用于图的生成树计数。
      
      看到给出图求生成树的这类问题就大概要往这方面想了。
      
      算法会根据图构造出一个特殊的基尔霍夫矩阵(A),接着根据矩阵树定理,用(A)计算出生成树个数。
      
      
      

    1.无向图的生成树计数

      
      对于给定的可含重边的连通无向图(G),求其生成树的个数。求法如下:
      
      定义度数矩阵(D):该矩阵仅在对角线上有值,(D_{i,i})表示(i)号点的度数。对于图中每一条无向边((u,v))(D_{u,u})++,(D_{v,v})++。
      
      定义邻接矩阵(C)(C_{i,j})表示(i)(j)的边数。对于图中每一条无向边((u,v))(C_{u,v})++,(C_{v,u})++。
      
      定义图(G)基尔霍夫矩阵(A=D-C)
      
    ​  矩阵树定理:将(A)去掉第(i)行和第(i)列((iin[1,n])),将它当做一个行列式求解,则(det(A))就是生成树个数。
      
        
      

    2.有向图的树形图计数

      
      对于有向图,不存在“生成树”的概念,但存在“树形图”的概念。有向图中,若选定一个点作为树根,能构造出一棵“树”(包含(n-1)条边)使得根能到达任意节点,则这是一棵外向树;若能构造出一棵“树”使得任意节点能到达根,则这是一棵内向树。
      
      定义度数矩阵(D):该矩阵仅在对角线上有值,(D_{i,i})表示(i)号点的度数。
      
               对于图中每一条有向边((u,v)),若构造外向树则(D_{v,v})++;若构造内向树则(D_{u,u})++。
      
      定义邻接矩阵(C)(C_{i,j})表示(i)(j)的边数。对于图中每一条有向边((u,v))(C_{u,v})++。
      
    ​  定义图(G)基尔霍夫矩阵(A=D-C)
      
      矩阵树定理:将(A)去掉第(i)行和第(i)列((iin[1,n])),将它当做一个行列式求解,则(det(A))就是以(i)为根的外向/内向树形图个数。很多时候我们会发现(A)的对角线上某数为(A_{i,i}=0),删去第(i)行和第(i)列可以干掉0。只有这样行列式才不等于0,其实也就是说只能从(i)出发有解了。
      
       
      

    3.细节

      
      求行列式的方法是:将行列式通过行列式初等变换消成上三角。此时对角线乘积即为行列式的值。
      
      注意,矩阵树定理这一套算法会考虑如何把所有的参与点构建成生成树,所以编号不要跳跃。如果说有障碍之类的元素,千万不要在矩阵中给它留一行一列,因为这一行一列都一定是0,算法会尝试将“障碍”构建进生成树,最后只能得到无解。
      
      
      
      
      

    一些例题

      

    BZOJ 4031

      
      传送门在此
      
      这是一道无向图生成树计数的裸题,直接上基础算法即可。
      
      这里就要注意障碍的处理了。我们应该对非障碍格子重新标号使得它们的编号连续,保证算法正常进行。
      
      这题的模数真的恶心,不是质数,没法用逆元。所以使用类辗转相除法来消元,每行的消元从(O(n))变成(O(nlg ))
      

    #include <cstdio>
    using namespace std;
    const int N=10,MOD=1e9;
    int n,m,id[N][N],idcnt,a[N*N][N*N];
    char map[N][N];
    inline void swap(int &x,int &y){x^=y^=x^=y;}
    inline int plus(int x,int y){return (x+y)%MOD;}
    inline int mul(int x,int y){return 1LL*x*y%MOD;}
    inline bool ok(int x,int y){return 1<=x&&x<=n&&1<=y&&y<=m&&map[x][y]=='.';}
    void addEdge(int u,int v){
        a[u][u]++; a[v][v]++;
        a[u][v]--; a[v][u]--;
    }
    int solve(){
        if(idcnt==1) return 1;
        int all=idcnt-1,res=1;
        for(int i=1;i<=all;i++)
            for(int j=i+1;j<=all;j++)
                while(a[j][i]){
                    int t=a[i][i]/a[j][i];
                    for(int k=i;k<=all;k++){
                        int q=plus(a[i][k],-mul(a[j][k],t));
                        a[i][k]=a[j][k];
                        a[j][k]=q;
                    }
                    res=-res;
                }
        for(int i=1;i<=all;i++) res=mul(res,a[i][i]);
        return (res+MOD)%MOD;
    }
    int main(){
        scanf("%d%d",&n,&m);
        for(int i=1;i<=n;i++) scanf("%s",map[i]+1);
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
                if(map[i][j]=='.') id[i][j]=++idcnt;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=m;j++)
                if(ok(i,j)){
                    if(ok(i+1,j))
                        addEdge(id[i][j],id[i+1][j]);
                    if(ok(i,j+1))
                        addEdge(id[i][j],id[i][j+1]);
                }
        printf("%d
    ",solve());
        return 0;
    }
    

      
      
      
      
      

    BZOJ 4894

      
      传送门
      
      这是一道有向图树形图计数。要求以1号点为根的外向树形图个数。
      
      按照上述做法直接写即可。删去(A)的第1行第1列,因为1号点没有入边,若不删第一行第一列行列式值为0,无法计算。
      

    #include <cstdio>
    using namespace std;
    const int N=305,MOD=1e9+7;
    int n,a[N][N];
    char str[N];
    inline int mul(int x,int y){return 1LL*x*y%MOD;}
    inline int plus(int x,int y){return (x+y)%MOD;}
    inline void swap(int &x,int &y){x^=y^=x^=y;}
    int ksm(int x,int y){
        int res=1;
        for(;y;x=mul(x,x),y>>=1)
            if(y&1) res=mul(res,x);
        return res;
    }
    int gaussian(){
        int res=1,size=n-1;
        for(int i=1;i<size;i++){
            if(!a[i][i]){
                int l;
                for(l=i+1;l<=size;l++)
                    if(a[l][i]) break;
                if(l<=size&&a[l][i]){
                    for(int j=i;j<=size;j++) swap(a[l][j],a[i][j]);
                    res=-res;
                }
                else return 0;
            }
            for(int j=i+1;j<=size;j++){
                int t=mul(a[j][i],ksm(a[i][i],MOD-2));
                for(int k=i;k<=size;k++)         
                    a[j][k]=plus(a[j][k],-mul(a[i][k],t));
            }
        }
        for(int i=1;i<=size;i++) res=mul(res,a[i][i]);
        return plus(res,MOD);
    }
    int main(){
        scanf("%d",&n);
        for(int i=1;i<=n;i++){
            scanf("%s",str+1);
            for(int j=1;j<=n;j++)
                if(str[j]=='1') 
                    a[j][j]++,a[i][j]--;
        }
        for(int i=1;i<n;i++)
            for(int j=1;j<n;j++) a[i][j]=a[i+1][j+1];
        printf("%d
    ",gaussian());
        return 0;
    }
    

      
      
      
      
      

    BZOJ 4596

      
      传送门
      
    ​  我开始不会做了。
      
    ​  我们发现如果将所有公司提供的边都加进图中,然后求生成树个数,是无法限制“每个公司至少要建一条”这个条件的。有的生成树可能只有一部分公司参与,比如说某一种生成树只含有(S)集合的公司。
      
      如果仅加入(S)集合所含公司的边,我们发现这些答案也会被统计到。既然有重复统计,可以考虑消除吗?
      
    ​  于是容斥的思想就体现出来了!
      
      我们枚举加入哪一些公司,分别求生成树个数。记参与公司集合为(S)时生成树个数为(f(S)),记有(x)个公司参与时的生成树总方案为(中有个A_x=sum_{S中有x个}f(S)),则

    [Ans=A_n-A_{n-1}+A_{n-2}-A_{n-3}...... ]

    ​  复杂度为(O(2^nn^3)),其实是可以跑的过的。
      

    #include <cstdio>
    #include <vector>
    #define mp make_pair
    #define pb push_back
    using namespace std;
    typedef pair<int,int> pii;
    const int N=18,MOD=1e9+7;
    int n;
    int a[N][N];
    vector<pii> l[N];
    vector<int> b[N];
    inline bool in(int i,int j){return (i>>(j-1))&1;}
    inline int plus(int x,int y){return (x+y)%MOD;}
    inline int mul(int x,int y){return 1LL*x*y%MOD;}
    inline void swap(int &x,int &y){x^=y^=x^=y;}
    void clear_mat(){
    	for(int i=1;i<=n;i++)
    		for(int j=1;j<=n;j++) a[i][j]=0;
    }
    void add_company(int id){
    	for(int i=0,sz=l[id].size();i<sz;i++){
    		int u=l[id][i].first,v=l[id][i].second;
    		a[u][u]++; a[v][v]++;
    		a[u][v]--; a[v][u]--;
    	}
    }
    int ksm(int x,int y){
    	int res=1;
    	for(;y;x=mul(x,x),y>>=1)
    		if(y&1) res=mul(res,x);
    	return res;
    }
    int gaussian(){
    	int s=n-1,res=1;
    	for(int i=1;i<=s;i++){
    		if(!a[i][i]){
    			int t;
    			for(t=i+1;t<=s&&!a[t][i];t++);
    			if(t>s) return 0;
    			for(int j=i;j<=s;j++) swap(a[i][j],a[t][j]);
    			res=-res;
    		}
    		int inv=ksm(a[i][i],MOD-2);
    		for(int j=i+1;j<=s;j++){
    			int t=mul(a[j][i],inv);
    			for(int k=i;k<=s;k++)
    				a[j][k]=plus(a[j][k],-mul(a[i][k],t));
    		}
    	}
    	for(int i=1;i<=s;i++) res=mul(res,a[i][i]);
    	return res;
    }
    int main(){
    	scanf("%d",&n);
    	for(int i=1,m;i<n;i++){
    		scanf("%d",&m);
    		for(int j=1,u,v;j<=m;j++){
    			scanf("%d%d",&u,&v);
    			l[i].pb(mp(u,v));	
    		}
    	}
    	int all=1<<(n-1);
    	for(int i=1;i<all;i++){
    		int cnt=0;
    		for(int j=1;j<=n;j++) 
    			cnt+=in(i,j);
    		b[cnt].pb(i);	
    	}
    	int ans=0;
    	for(int i=n-1,r=1;i>=1;i--,r=-r)
    		for(int j=0,sz=b[i].size();j<sz;j++){
    			clear_mat();
    			int st=b[i][j];
    			for(int c=1;c<n;c++)
    				if(in(st,c)) add_company(c);
    			ans=plus(ans,gaussian()*r);				
    		}
    	ans=plus(ans,MOD);
    	printf("%d
    ",ans);
    	return 0;
    }
    

      
      
      

    总结

      
      矩阵树定理本身还是挺简单的,但愿自己不要忘得太快......
      
      但是要灵活运用(废话)
      
    ​  如果要深入透彻,我还是得研究一下矩阵树定理的证明。不过就当一个大坑先留着吧。

  • 相关阅读:
    MyBatis的动态SQL语句这么厉害的!
    连接数据库,使用c3p0技术连接MySQL数据库
    Servlet 常见的乱码解决方案
    超级签具体实现
    Xcode报错You don’t have permission.
    SpringBoot+Mybatis整合实例
    恢复mysql数据库误删数据
    日期(date)运用座谈会
    程序猿日记--学习怎样学习
    服务器数据库密码忘记
  • 原文地址:https://www.cnblogs.com/RogerDTZ/p/9226071.html
Copyright © 2011-2022 走看看