矩阵树定理
矩阵树定理
前言：欢迎来到线性代数最高城——Matrix-Tree定理。
写在矩阵树定理之前
行列式
定义一个\(N\times N\)的矩阵\(A\)。
定义一个矩阵的行列式：
其中\(P\)为一个\(1\sim n\)的排列，\(k\)为\(P\)的逆序对数。
如果直接做，容易发现时间复杂度为\(O(N!\cdot N)\)。
考虑一些特殊矩阵，比如：
当\(A=I\)，\(I\)为单位矩阵时，总共的排列数为\(N!\)，当且仅当\(P=\{1,2,3,\cdots n\}\)时才会对答案产生贡献。P.S.其他的情况中间会乘上\(0\)导致没有贡献。
此时，\(k=0\)，所以\(|I|=(-1)^0=1\)。
考虑使用高斯消元法将\(A\)转换为单位矩阵\(I\)。对于高斯消元的每次操作，我们都可以将其看做对原矩阵做矩阵操作，我们将操作到过来就能得到\(A\)的行列式。
分别考虑高斯消元操作总对答案贡献的影响，设变换后的矩阵为\(A^{'}\)。
- 交换\(A\)中的两行。
容易发现加法乘法都有交换律，所以对于一个排列\(P\)对应的\(A_{P_i}\)，一定能在操作后的矩阵中找到相同的的数组。
更具体的，在数组\(A_{P_i}\)交换的两行中，设它们为\(A_{P_i},A_{P_j}\)。
交换后的数组\(A_{P_i}^{'}\)，我们同时交换\(P_i\)，所以有\(A^{'}_{P^{'}_i}=A_{P_i},A^{'}_{P^{'}_j}=A_{P_j}\)
显然，这是双射的，意思是这是一一对应的。
此时，每一个排列的逆序对数都乘上/除以\(-1\)，此时\(|A^{'}|=|A|\)。
- 将其中一行乘上\(K\)。
此时每一个排列都有且仅有一个元素变为原来的\(K\)倍。
我们将\(K\)提出来，有\(|A^{'}|=K\cdot|A|\)。
- 如果有两行完全相同，那么\(|A|=0\)。
运用性质1，我们有\(|A|=-|A^{'}|\)，由于这两行完全相同，我们又有\(A=A^{'}\)。
所以\(|A|=-|A|\)，所以\(|A|=0\)。
通过性质2推广，如果一行正好是另一行的\(K\)倍，那么有\(|A|=|A^{'}|\)
- 将其中一行加上另一行的\(K\)倍。
我们拆式子，发现我们相当于加上的一个新矩阵的行列式，这个新矩阵有一行正好是另一行的\(K\)倍。
所以新矩阵的贡献为\(0\)，所以\(|A|=|A^{'}|\)。
所以我们做一次高斯消元，同时维护行列式的值。容易发现，\(|A|\)的操作只有取反或\(\times K\)，将\(\times K\)变为乘上\(K\)的逆元，这个操作具有交换律，所以我们将初值设为\(1\)就行了。
矩阵树定理
给一张有向图\(G=(V,E)\)，定义入度/出度矩阵\(D_{in/out}\)，在\(i=j\)时为该点的入度/出度，其余全为\(0\)。
再定义邻接矩阵\(A\)，在\((i,j)\in E\)是时为\(1\)，其余为\(0\)。
矩阵树定理：以\(i\)为根的外向树个数为\(D_{in}-A\)删去第\(i\)行和第\(i\)列的行列式的值。以\(i\)为根的内向树个数为\(D_{out}-A\)删去第\(i\)行和第\(i\)列的行列式的值。
当\(G\)有边权时，我们更改\(A\)的定义，我们将\((i,j)\in E\)时权值为边权，其余为\(0\)，同时，设每个点的入度/出度为入边/出边的权值和。
这时每棵生成树的权值为边权的乘积。
证明不会。。。不过这里有可能详细的证明。
使用高斯消元法求解的时间复杂度为\(O(N^3)\)。
应用
Solution
维护以一为根的外向生成树的权值，删去第一行第一列。
使用高斯消元法维护矩阵的行列式的值，时间复杂度为\(O(N^3)\)。
代码如下：
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N=310,mod=1e9+7;
ll power(ll a,ll b){
ll ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
int n,m,op;
ll A[N][N];
ll solve(){
int ans=1,opt=1;
for(int i=2,pos;i<=n;i++){
pos=i;
for(int j=i+1;j<=n;j++)
if(A[j][i]>A[pos][i])pos=j;
if(i!=pos)opt=-opt,swap(A[i],A[pos]);
if(!A[i][i])return 0;
for(int j=2;j<=n;j++){
if(j==i)continue;
ll tmp=A[j][i]*power(A[i][i],mod-2)%mod;
for(int k=2;k<=n;k++)
A[j][k]=(A[j][k]-tmp*A[i][k]%mod+mod)%mod;
}
ans=ans*A[i][i]%mod;
ll tmp=power(A[i][i],mod-2);
for(int j=2;j<=n;j++)A[i][j]=A[i][j]*tmp%mod;
}
return (ans*opt%mod+mod)%mod;
}
int main(){
scanf("%d%d%d",&n,&m,&op);
for(int i=1,x,y,z;i<=m;i++){
scanf("%d%d%d",&x,&y,&z);
if(op==0){
A[x][x]=(A[x][x]+z)%mod,A[y][y]=(A[y][y]+z)%mod;
A[x][y]=(A[x][y]-z+mod)%mod,A[y][x]=(A[y][x]-z+mod)%mod;
}else if(op==1){
A[y][y]=(A[y][y]+z)%mod;
A[x][y]=(A[x][y]-z+mod)%mod;
}
}
printf("%lld\n",solve());
return 0;
}
题意简述：每条边概率连通，求最求概率最小生成树。
Solution
由上面的例题我们容易发现，矩阵树定理求的是对于图上所有的生成树\(T\)，设\(w(e)\)为\(e\)的边权，求解：
通过概率的加法/乘法原理，我们将所有的生成树的概率加起来，一颗生成树的概率为\(\prod_{e\in T}P_e\prod_{e\notin T}(1-P_e)\)，所以答案就是：
推推式子，将式子化为\(\sum_T\prod_{e\in T}w(e)\)的样子。
容易发现\(\prod_{e\in G}(1-P_e)\)是定值，将它移到外面。
我们将边权设为\(\frac{P_e}{(1-Pe)}\)，跑矩阵树定理即可。
代码如下：
#include<bits/stdc++.h>
using namespace std;
const int N=70;
const double eps=1e-7;
double A[N][N],ans=1.0;
int n;
double solve(){
double res=1.0;
int opt=1;
for(int i=1;i<=n-1;i++){
int pos=i;
for(int j=i+1;j<=n-1;j++)
if(A[pos][i]<A[j][i])pos=j;
if(pos!=i)swap(A[i],A[pos]),opt=-opt;
if(fabs(A[i][i])<eps)return 0.0;
for(int j=i+1;j<=n-1;j++){
double tmp=A[j][i]/A[i][i];
for(int k=i;k<=n-1;k++)
A[j][k]-=tmp*A[i][k];
}
res=res*A[i][i];
for(int j=i;j<=n-1;j++)A[i][j]=A[i][j]/A[i][i];
}
return res*opt;
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
scanf("%lf",&A[i][j]);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++){
if(fabs(A[i][j])<eps)A[i][j]=eps;
else if(fabs(1.0-A[i][j])<eps)A[i][j]=1.0-eps;
if(j>i)ans*=(1-A[i][j]);
A[i][j]=A[i][j]/(1-A[i][j]);
}
for(int i=1;i<=n;i++){
A[i][i]=0;
for(int j=1;j<=n;j++)if(i!=j){
A[i][i]+=A[i][j],A[i][j]=-A[i][j];
}
}
ans*=solve();
printf("%.10lf\n",ans);
return 0;
}
给定一张图，边有颜色，求正好\(N-1\)个颜色的生成树。
Solution
容易发现，我们强制让\(N-1\)个颜色出现不太好做，但是我们至多\(K\)个颜色是很好做的。
我们\(K\)个颜色所属的边权加在一张图中，跑矩阵树定理就做完了。
于是，我们要让一个“至多”的条件变为严格“等于”的条件。
由于\(N\le 17\)，所以我们直接\(O(2^N)\)容斥即可。
时间复杂度为\(O(N^32^N)\)。
代码如下：
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N=30,mod=1e9+7;
ll power(ll a,ll b){
ll ans=1;
for(;b;b>>=1){
if(b&1)ans=ans*a%mod;
a=a*a%mod;
}
return ans;
}
ll A[N][N],ans;
int n;
vector<pair<ll,ll>>g[N];
ll solve(){
ll res=1,opt=1;
for(int i=1;i<=n;i++){
int pos=-1;
for(int j=i;j<=n;j++)
if(A[j][i]){pos=j;break;}
if(pos==-1)return 0;
if(i!=pos)opt=-opt,swap(A[i],A[pos]);
for(int j=1;j<=n;j++){
if(i==j)continue;
ll tmp=A[j][i]*power(A[i][i],mod-2)%mod;
for(int k=1;k<=n;k++)
A[j][k]=(A[j][k]-tmp*A[i][k]%mod+mod)%mod;
}
res=res*A[i][i]%mod;
ll tmp=power(A[i][i],mod-2);
for(int j=1;j<=n;j++)A[i][j]=A[i][j]*tmp%mod;
}
return (res*opt%mod+mod)%mod;
}
int main(){
scanf("%d",&n),n--;
for(int i=1,K;i<=n;i++){
scanf("%d",&K);
while(K--){
int x,y;scanf("%d%d",&x,&y);
g[i].push_back({x,y});
}
}
for(int i=1;i<1<<n;i++){
memset(A,0,sizeof(A));
int cnt=0;
for(int j=1;j<=n;j++){
if(!(i>>(j-1)&1))continue;
cnt++;
for(int k=0;k<g[j].size();k++){
int x=g[j][k].first,y=g[j][k].second;
A[x][y]++,A[y][x]++,A[x][x]--,A[y][y]--;
}
}
ans=((ans+((cnt&1)?-1:1)*solve())%mod+mod)%mod;
}
printf("%lld\n",ans);
return 0;
}
浙公网安备 33010602011771号