【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\)列（\(i\in[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\)列（\(i\in[1,n]\)），将它当做一个行列式求解，则\(\det(A)\)就是以\(i\)为根的外向/内向树形图个数。很多时候我们会发现\(A\)的对角线上某数为\(A_{i,i}=0\)，删去第\(i\)行和第\(i\)列可以干掉0。只有这样行列式才不等于0，其实也就是说只能从\(i\)出发有解了。
3.细节
求行列式的方法是：将行列式通过行列式初等变换消成上三角。此时对角线乘积即为行列式的值。
注意，矩阵树定理这一套算法会考虑如何把所有的参与点构建成生成树，所以编号不要跳跃。如果说有障碍之类的元素，千万不要在矩阵中给它留一行一列，因为这一行一列都一定是0，算法会尝试将“障碍”构建进生成树，最后只能得到无解。
一些例题
BZOJ 4031
传送门在此
这是一道无向图生成树计数的裸题，直接上基础算法即可。
这里就要注意障碍的处理了。我们应该对非障碍格子重新标号使得它们的编号连续，保证算法正常进行。
这题的模数真的恶心，不是质数，没法用逆元。所以使用类辗转相除法来消元，每行的消元从\(O(n)\)变成\(O(n\lg )\)。
#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\n",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\n",gaussian());
return 0;
}
BZOJ 4596
传送门
我开始不会做了。
我们发现如果将所有公司提供的边都加进图中，然后求生成树个数，是无法限制“每个公司至少要建一条”这个条件的。有的生成树可能只有一部分公司参与，比如说某一种生成树只含有\(S\)集合的公司。
如果仅加入\(S\)集合所含公司的边，我们发现这些答案也会被统计到。既然有重复统计，可以考虑消除吗？
于是容斥的思想就体现出来了！
我们枚举加入哪一些公司，分别求生成树个数。记参与公司集合为\(S\)时生成树个数为\(f(S)\)，记有\(x\)个公司参与时的生成树总方案为\(中有个A_x=\sum_{S中有x个}f(S)\)，则
复杂度为\(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\n",ans);
return 0;
}
总结
矩阵树定理本身还是挺简单的，但愿自己不要忘得太快......
但是要灵活运用（废话）。
如果要深入透彻，我还是得研究一下矩阵树定理的证明。不过就当一个大坑先留着吧。
