引言
矩阵树定理是主要为了解决生成树计数的问题，生成树计数分为无向图计数以及有向图计数。除此之外，变元矩阵树定理还能够求解生成树边权积的和。
本文会先给出矩阵树定理的基本和扩展内容，然后给出关于定理的详细证明，最后是相关的习题及讲解。
定理内容
生成树计数
无向图
对于一般无向图而言，将图中的自环的边去除以后，计算出满足如下条件的基尔霍夫(Kirchhoff)矩阵：
L
i
,
j
=
{
i
点
度
数
,
if
i
=
j
i
与
j
连
边
数
的
相
反
数
,
if
i
≠
j
\mathbf{L_{i,j} = \begin{cases} i点度数, & \text{if }i=j \\ i与j连边数的相反数, & \text{if }i\ne j \end{cases}}
Li,j={i点度数,i与j连边数的相反数,if i=jif i=j
那么生成树的个数就等于矩阵
L
\mathbf{L}
L的任意一个代数余子式，一般去掉矩阵
L
\mathbf{L}
L的最后一行和最后一列后作为答案。
有向图
对于有向图而言，定理需要做一些改变，有向图的树形图分为两种，一是外向树(即从祖先节点连向其子孙)，二是内向树(即从子孙连向祖先节点)，定理改变如下：
外向树：
L
i
,
j
′
=
{
i
点
入
度
,
if
i
=
j
i
点
连
向
j
点
的
边
数
的
相
反
数
,
if
i
≠
j
\mathbf{L'_{i,j} = \begin{cases} i点入度, & \text{if }i=j \\ i点连向j点的边数的相反数, & \text{if }i\ne j \end{cases}}
Li,j′={i点入度,i点连向j点的边数的相反数,if i=jif i=j
内向树：
L
i
,
j
′
′
=
{
i
点
出
度
,
if
i
=
j
i
点
连
向
j
点
的
边
数
的
相
反
数
,
if
i
≠
j
\mathbf{L''_{i,j} = \begin{cases} i点出度, & \text{if }i=j \\ i点连向j点的边数的相反数, & \text{if }i\ne j \end{cases}}
Li,j′′={i点出度,i点连向j点的边数的相反数,if i=jif i=j
以i为根节点的树形图个数就是对应矩阵去除第i行和第i列后的行列式。
定理扩展
除了生成树计数以外，变元的矩阵树定理还可以求解生成树边权乘积的和。
变元矩阵树定理需要对基尔霍夫矩阵作一些改变，改变后的基尔霍夫矩阵如下：(与上文的三个公式对应，矩阵用ℒ符号表示)
L
i
,
j
=
{
与
i
相
连
边
权
和
,
if
i
=
j
i
与
j
连
边
权
和
的
相
反
数
,
if
i
≠
j
\mathbf{ℒ_{i,j} = \begin{cases} 与i相连边权和, & \text{if }i=j \\ i与j连边权和的相反数, & \text{if }i\ne j \end{cases}}
Li,j={与i相连边权和,i与j连边权和的相反数,if i=jif i=j
L i , j ′ = { 以 i 点 结 束 边 权 和 , if i = j i 点 连 向 j 点 的 边 权 和 的 相 反 数 , if i ≠ j \mathbf{ℒ'_{i,j} = \begin{cases} 以i点结束边权和, & \text{if }i=j \\ i点连向j点的边权和的相反数, & \text{if }i\ne j \end{cases}} Li,j′={以i点结束边权和,i点连向j点的边权和的相反数,if i=jif i=j
L
i
,
j
′
′
=
{
从
i
点
出
发
边
权
和
,
if
i
=
j
i
点
连
向
j
点
的
边
权
和
的
相
反
数
,
if
i
≠
j
\mathbf{ℒ''_{i,j} = \begin{cases} 从i点出发边权和, & \text{if }i=j \\ i点连向j点的边权和的相反数, & \text{if }i\ne j \end{cases}}
Li,j′′={从i点出发边权和,i点连向j点的边权和的相反数,if i=jif i=j
计算行列式的方法与原定理对应相同，对于图G而言，计算出来的结果用公式表示出来如下：
∑
T
r
e
e
∈
G
∏
e
d
g
e
∈
T
r
e
e
v
a
l
[
e
d
g
e
]
\mathbf{\sum_{Tree\in G}\prod_{edge\in Tree}val\lbrack edge\rbrack}
Tree∈G∑edge∈Tree∏val[edge]
定理证明
待补充。。。
定理应用
例题一
题目来源：洛谷【模板】Matrix-Tree 定理
题解：直接套用公式就行了，不过需要注意下行列式求解的过程，我是直接将行列式的第一行置为1，第一列置为0(除了
L
1
,
1
\mathbf{L_{1,1}}
L1,1)，这等价于代数余子式。
代码：
#include <bits/stdc++.h>
#define FOR(i,a,b) for(register int i=(a);i<(b);++i)
#define ROF(i,a,b) for(register int i=(a);i>=(b);--i)
#define pi pair<int,int>
#define mk(a,b) make_pair(a,b)
#define mygc(c) (c)=getchar()
#define mypc(c) putchar(c)
#define fi first
#define se second
using namespace std;
typedef long long ll;
typedef double db;
const int maxn = 305;
const int maxm = 20000;
const int inf = 2147483647;
typedef long long ll;
const double eps = 1e-9;
const long long INF = 9223372036854775807ll;
ll qpow(ll a,ll b,ll c){ll ans=1;while(b){if(b&1)ans=ans*a%c;a=a*a%c;b>>=1;}return ans;}
inline void rd(int *x){int k,m=0;*x=0;for(;;){mygc(k);if(k=='-'){m=1;break;}if('0'<=k&&k<='9'){*x=k-'0';break;}}for(;;){mygc(k);if(k<'0'||k>'9')break;*x=(*x)*10+k-'0';}if(m)(*x)=-(*x);}
inline void rd(ll *x){int k,m=0;*x=0;for(;;){mygc(k);if(k=='-'){m=1;break;}if('0'<=k&&k<='9'){*x=k-'0';break;}}for(;;){mygc(k);if(k<'0'||k>'9')break;*x=(*x)*10+k-'0';}if(m)(*x)=-(*x);}
inline void rd(db *x){scanf("%lf",x);}
inline int rd(char c[]){int i,s=0;for(;;){mygc(i);if(i!=' '&&i!='\n'&&i!='\r'&&i!='\t'&&i!=EOF) break;}c[s++]=i;for(;;){mygc(i);if(i==' '||i=='\n'||i=='\r'||i=='\t'||i==EOF) break;c[s++]=i;}c[s]='\0';return s;}
inline void rd(int a[],int n){FOR(i,0,n)rd(&a[i]);}
inline void rd(ll a[],int n){FOR(i,0,n)rd(&a[i]);}
template <class T, class S> inline void rd(T *x, S *y){rd(x);rd(y);}
template <class T, class S, class U> inline void rd(T *x, S *y, U *z){rd(x);rd(y);rd(z);}
template <class T, class S, class U, class V> inline void rd(T *x, S *y, U *z, V *w){rd(x);rd(y);rd(z);rd(w);}
inline void wr(int x){if(x < 10) putchar('0' + x); else wr(x / 10), wr(x % 10);}
inline void wr(int x, char c){int s=0,m=0;char f[10];if(x<0)m=1,x=-x;while(x)f[s++]=x%10,x/=10;if(!s)f[s++]=0;if(m)mypc('-');while(s--)mypc(f[s]+'0');mypc(c);}
inline void wr(ll x, char c){int s=0,m=0;char f[20];if(x<0)m=1,x=-x;while(x)f[s++]=x%10,x/=10;if(!s)f[s++]=0;if(m)mypc('-');while(s--)mypc(f[s]+'0');mypc(c);}
inline void wr(db x, char c){printf("%.15f",x);mypc(c);}
inline void wr(const char c[]){int i;for(i=0;c[i]!='\0';i++)mypc(c[i]);}
inline void wr(const char x[], char c){int i;for(i=0;x[i]!='\0';i++)mypc(x[i]);mypc(c);}
template<class T> inline void wrn(T x){wr(x,'\n');}
template<class T, class S> inline void wrn(T x, S y){wr(x,' ');wr(y,'\n');}
template<class T, class S, class U> inline void wrn(T x, S y, U z){wr(x,' ');wr(y,' ');wr(z,'\n');}
template<class T, class S, class U,class H> inline void wrn(T x, S y, U z,H h){wr(x,' ');wr(y,' ');wr(z,' ');wr(h,'\n');}
template<class T> inline void wra(T x[], int n){int i;if(!n){mypc('\n');return;}FOR(i,0,n-1)wr(x[i],' ');wr(x[n-1],'\n');}
const ll mod=1e9+7;
struct Matrix{
ll a[maxn][maxn],n,m;
Matrix(int n=0,int m=0):n(n),m(m){}
void clear(){
memset(a,0,sizeof(a));
}
void I(){
FOR(i,0,n)
FOR(j,0,m)if(i==j)a[i][j]=1;else a[i][j]=0;
}
void init(ll b[],int nn,int mm){
n=nn,m=mm;
FOR(i,0,n)
FOR(j,0,m)a[i][j]=b[i*mm+j];
}
ll det(){//求解取模情况下的det
int maxrow=0,fg=1;
ll temp;
FOR(i,0,n){
maxrow=i;
FOR(j,i,n)if(a[j][i]>a[i][i])maxrow=j;
if(maxrow!=i)swaprow(i,maxrow,i),fg*=-1;
ll reva=qpow(a[i][i],mod-2,mod);
FOR(j,i+1,n){
temp=a[j][i]*reva%mod;
FOR(k,i,m){
a[j][k]=(a[j][k]-a[i][k]*temp%mod)%mod;
a[j][k]=(a[j][k]%mod+mod)%mod;
}
}
}
ll ans=1;
FOR(i,0,n)ans=ans*a[i][i]%mod;
ans*=fg;
ans=(ans%mod+mod)%mod;
return ans;
}
void swaprow(int r1,int r2,int beg){
FOR(i,beg,m)swap(a[r1][i],a[r2][i]);
return;
}
void print(){
FOR(i,0,n){
FOR(j,0,m){
printf("%lld ",a[i][j]);
}
printf("\n");
}
printf("\n");
}
void printrow(int r){
FOR(i,0,m)printf("%lld ",a[r][i]);
printf("\n");
}
void printcol(int c){
FOR(i,0,n)printf("%lld ",a[i][c]);
printf("\n");
}
}c;
int main(){
int n,m,t;
rd(&n,&m,&t);
c.n=c.m=n;
while(m--){
int u,v,w;
rd(&u,&v,&w);
u--,v--;
if(!t){
c.a[u][v]=(c.a[u][v]-w+mod)%mod;
c.a[v][u]=(c.a[v][u]-w+mod)%mod;
c.a[u][u]=(c.a[u][u]+w)%mod,c.a[v][v]=(c.a[v][v]+w)%mod;
}else{
c.a[u][v]=(c.a[u][v]-w+mod)%mod;
c.a[v][v]=(c.a[v][v]+w)%mod;
}
}
// c.print();
FOR(i,0,n)c.a[0][i]=1,c.a[i][0]=0;
c.a[0][0]=1;
wrn(c.det());
// c.print();
}
例题二
题目来源：洛谷[SDOI2014]重建
题面：
题解：本题主要考察对变元矩阵树定理公式的熟悉程度。
先通过基本的乘法原理得出答案，对于每一个生成树而言都一定能够保证所有城市连通，所以形成生成树的概率是：
∑
T
r
e
e
∈
G
∏
e
d
g
e
∈
T
r
e
e
p
[
e
d
g
e
]
\mathbf{\sum_{Tree\in G}{\prod_{edge\in Tree}{p[ edge]}}}
Tree∈G∑edge∈Tree∏p[edge]
但是保证生成树的存在还不够，还要保证生成树以外的边都断开，根据乘法原理只需要乘上其他边断开的概率即可：
∏
e
d
g
e
∉
T
r
e
e
(
1
−
p
[
e
d
g
e
]
)
⋅
∑
T
r
e
e
∈
G
∏
e
d
g
e
∈
T
r
e
e
p
[
e
d
g
e
]
\mathbf{\prod_{edge\notin Tree}(1-p[edge]) \cdot \sum_{Tree\in G}\prod_{edge\in Tree}p[edge]}
edge∈/Tree∏(1−p[edge])⋅Tree∈G∑edge∈Tree∏p[edge]
下面我们做一些变形，对整体乘一个
∏
e
d
g
e
∈
T
r
e
e
(
1
−
p
[
e
d
g
e
]
)
\mathbf{\prod_{edge\in Tree}(1-p[edge]) }
∏edge∈Tree(1−p[edge])，再除回去就得到：
原
式
=
∏
e
d
g
e
∈
T
r
e
e
(
1
−
p
[
e
d
g
e
]
)
⋅
∏
e
d
g
e
∉
T
r
e
e
(
1
−
p
[
e
d
g
e
]
)
⋅
∑
T
r
e
e
∈
G
∏
e
d
g
e
∈
T
r
e
e
p
[
e
d
g
e
]
∏
e
d
g
e
∈
T
r
e
e
(
1
−
p
[
e
d
g
e
]
)
=
∏
e
d
g
e
∈
G
(
1
−
p
[
e
d
g
e
]
)
⋅
∑
T
r
e
e
∈
G
∏
e
d
g
e
∈
T
r
e
e
p
[
e
d
g
e
]
∏
e
d
g
e
∈
T
r
e
e
(
1
−
p
[
e
d
g
e
]
)
=
∏
e
d
g
e
∈
G
(
1
−
p
[
e
d
g
e
]
)
⋅
∑
T
r
e
e
∈
G
∏
e
d
g
e
∈
T
r
e
e
p
[
e
d
g
e
]
1
−
p
[
e
d
g
e
]
\mathbf{\begin{aligned} 原式 & = \prod_{edge\in Tree}(1-p[edge])\cdot \prod_{edge\notin Tree}(1-p[edge]) \cdot \frac{\sum_{Tree\in G}\prod_{edge\in Tree}p[edge]}{\prod_{edge\in Tree}(1-p[edge])} \\ & = \prod_{edge\in G}(1-p[edge]) \cdot \sum_{Tree\in G}{\frac{\prod_{edge\in Tree}p[edge]}{\prod_{edge\in Tree}(1-p[edge])}} \\ & = \prod_{edge\in G}(1-p[edge]) \cdot \sum_{Tree\in G}{\prod_{edge\in Tree}{\frac{p[edge]}{1-p[edge]}}} \\ \end{aligned} }
原式=edge∈Tree∏(1−p[edge])⋅edge∈/Tree∏(1−p[edge])⋅∏edge∈Tree(1−p[edge])∑Tree∈G∏edge∈Treep[edge]=edge∈G∏(1−p[edge])⋅Tree∈G∑∏edge∈Tree(1−p[edge])∏edge∈Treep[edge]=edge∈G∏(1−p[edge])⋅Tree∈G∑edge∈Tree∏1−p[edge]p[edge]
然后很容易观察到变形后的式子的右边项可以拿变元矩阵树定理来算，相当于只是将原定理中的
v
a
l
[
e
d
g
e
]
\mathbf{val[edge]}
val[edge]替换成了
p
[
e
d
g
e
]
1
−
p
[
e
d
g
e
]
\mathbf{\frac{p[edge]}{1-p[edge]}}
1−p[edge]p[edge]；而左边项在读入边的时候就可以提前算出来，最后两部分想乘即可。
最后需要注意的是 p [ e d g e ] = 1 \mathbf{p[edge]=1} p[edge]=1的情况，这种情况下分母为0，一个小技巧是加上一个很小的数即可。
代码：
#include <bits/stdc++.h>
#define FOR(i,a,b) for(register int i=(a);i<(b);++i)
#define ROF(i,a,b) for(register int i=(a);i>=(b);--i)
#define pi pair<int,int>
#define mk(a,b) make_pair(a,b)
#define mygc(c) (c)=getchar()
#define mypc(c) putchar(c)
#define fi first
#define se second
using namespace std;
typedef long long ll;
typedef double db;
const int maxn = 60;
const int maxm = 20000;
const int inf = 2147483647;
typedef long long ll;
const double eps = 1e-7;
const long long INF = 9223372036854775807ll;
ll qpow(ll a,ll b,ll c){ll ans=1;while(b){if(b&1)ans=ans*a%c;a=a*a%c;b>>=1;}return ans;}
inline void rd(int *x){int k,m=0;*x=0;for(;;){mygc(k);if(k=='-'){m=1;break;}if('0'<=k&&k<='9'){*x=k-'0';break;}}for(;;){mygc(k);if(k<'0'||k>'9')break;*x=(*x)*10+k-'0';}if(m)(*x)=-(*x);}
inline void rd(ll *x){int k,m=0;*x=0;for(;;){mygc(k);if(k=='-'){m=1;break;}if('0'<=k&&k<='9'){*x=k-'0';break;}}for(;;){mygc(k);if(k<'0'||k>'9')break;*x=(*x)*10+k-'0';}if(m)(*x)=-(*x);}
inline void rd(db *x){scanf("%lf",x);}
inline int rd(char c[]){int i,s=0;for(;;){mygc(i);if(i!=' '&&i!='\n'&&i!='\r'&&i!='\t'&&i!=EOF) break;}c[s++]=i;for(;;){mygc(i);if(i==' '||i=='\n'||i=='\r'||i=='\t'||i==EOF) break;c[s++]=i;}c[s]='\0';return s;}
inline void rd(int a[],int n){FOR(i,0,n)rd(&a[i]);}
inline void rd(ll a[],int n){FOR(i,0,n)rd(&a[i]);}
template <class T, class S> inline void rd(T *x, S *y){rd(x);rd(y);}
template <class T, class S, class U> inline void rd(T *x, S *y, U *z){rd(x);rd(y);rd(z);}
template <class T, class S, class U, class V> inline void rd(T *x, S *y, U *z, V *w){rd(x);rd(y);rd(z);rd(w);}
inline void wr(int x){if(x < 10) putchar('0' + x); else wr(x / 10), wr(x % 10);}
inline void wr(int x, char c){int s=0,m=0;char f[10];if(x<0)m=1,x=-x;while(x)f[s++]=x%10,x/=10;if(!s)f[s++]=0;if(m)mypc('-');while(s--)mypc(f[s]+'0');mypc(c);}
inline void wr(ll x, char c){int s=0,m=0;char f[20];if(x<0)m=1,x=-x;while(x)f[s++]=x%10,x/=10;if(!s)f[s++]=0;if(m)mypc('-');while(s--)mypc(f[s]+'0');mypc(c);}
inline void wr(db x, char c){printf("%.15f",x);mypc(c);}
inline void wr(const char c[]){int i;for(i=0;c[i]!='\0';i++)mypc(c[i]);}
inline void wr(const char x[], char c){int i;for(i=0;x[i]!='\0';i++)mypc(x[i]);mypc(c);}
template<class T> inline void wrn(T x){wr(x,'\n');}
template<class T, class S> inline void wrn(T x, S y){wr(x,' ');wr(y,'\n');}
template<class T, class S, class U> inline void wrn(T x, S y, U z){wr(x,' ');wr(y,' ');wr(z,'\n');}
template<class T, class S, class U,class H> inline void wrn(T x, S y, U z,H h){wr(x,' ');wr(y,' ');wr(z,' ');wr(h,'\n');}
template<class T> inline void wra(T x[], int n){int i;if(!n){mypc('\n');return;}FOR(i,0,n-1)wr(x[i],' ');wr(x[n-1],'\n');}
const ll mod=1e9+7;
struct Matrix{
db a[maxn][maxn];
int n,m;
Matrix(int n=0,int m=0):n(n),m(m){}
void clear(){
memset(a,0,sizeof(a));
}
void I(){
FOR(i,0,n)
FOR(j,0,m)if(i==j)a[i][j]=1;else a[i][j]=0;
}
void init(db b[],int nn,int mm){
n=nn,m=mm;
FOR(i,0,n)
FOR(j,0,m)a[i][j]=b[i*mm+j];
}
int dcmp(db x){
if(fabs(x)<eps)return 0;
if(x<0)return -1;
return 1;
}
void print(){
FOR(i,0,n){
FOR(j,0,m){
printf("%.4f ",a[i][j]);
}
printf("\n");
}
printf("\n");
}
void printrow(int r){
FOR(i,0,m)printf("%.4f ",a[r][i]);
printf("\n");
}
void printcol(int c){
FOR(i,0,n)printf("%.4f ",a[i][c]);
printf("\n");
}
int Gauss(){
tra();
int rmx=get_nzero_row();
if(rmx==-1)return 1;//特殊情况,多解
bool ok=false;
FOR(i,0,m-1)if(dcmp(a[rmx][i])){
ok=true;
break;
}
if(!ok)return -1;//无解
int res=-1;
if(rmx<m-2){//多解
FOR(i,rmx+1,m-1)a[i][m-1]=1;
ROF(i,rmx,0){
FOR(j,i+1,m-1)a[i][m-1]-=a[i][j]*a[j][m-1];
a[i][m-1]/=a[i][i];
}
res=1;
}else if(rmx==m-2){//一解
res=0;
ROF(i,rmx,0){
FOR(j,i+1,m-1)a[i][m-1]-=a[i][j]*a[j][m-1];
a[i][m-1]/=a[i][i];
}
}
return res;
}
void tra(){//转化为上三角矩阵
int maxrow=0;
db temp;
FOR(i,0,n){
maxrow=i;
FOR(j,i,n)if(fabs(a[j][i])>fabs(a[i][i]))maxrow=j;
if(maxrow!=i)swaprow(i,maxrow,i);
FOR(j,i+1,n){
temp=a[j][i]/a[i][i];
FOR(k,i,m){
a[j][k]-=a[i][k]*temp;
}
}
}
}
db det(){
int maxrow=0,fg=1;
db temp;
FOR(i,0,n){
maxrow=i;
FOR(j,i,n)if(fabs(a[j][i])>fabs(a[i][i]))maxrow=j;
if(maxrow!=i)swaprow(i,maxrow,i),fg*=-1;
FOR(j,i+1,n){
temp=a[j][i]/a[i][i];
FOR(k,i,m){
a[j][k]-=a[i][k]*temp;
}
}
}
db ans=1;
FOR(i,0,n)ans=ans*a[i][i];
ans*=fg;
return ans;
}
int get_nzero_row(){//返回最后一个非零行
ROF(i,n-1,0){
FOR(j,0,m)if(dcmp(a[i][j]))return i;
}
return -1;
}
void swaprow(int r1,int r2,int beg){
FOR(i,beg,m)swap(a[r1][i],a[r2][i]);
return;
}
}c;
int main(){
int n;
rd(&n);
c.n=c.m=n;
db ans=1.0;
FOR(i,0,n){
FOR(j,0,n){
db w,e;
rd(&w);
if(j<=i)continue;
if(fabs(1-w)<eps)w+=eps;
e=w/(1-w);
ans*=1-w;
c.a[i][i]+=e;
c.a[j][j]+=e;
c.a[i][j]=c.a[j][i]=-e;
}
}
FOR(i,0,n)c.a[0][i]=1.0,c.a[i][0]=0;
c.a[0][0]=1;
ans*=c.det();
wrn(ans);
}
例题三
待补充。。。
2374
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
