矩阵树定理

· · 个人记录

矩阵树定理

前言：欢迎来到线性代数最高城——Matrix-Tree定理。

写在矩阵树定理之前

行列式

定义一个N\times N的矩阵A。

定义一个矩阵的行列式：

|A|=\sum_P(-1)^k\prod_{i=1}^nA_{i,P_i}

其中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^{'}。

  1. 交换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|。

  1. 将其中一行乘上K。

此时每一个排列都有且仅有一个元素变为原来的K倍。

我们将K提出来，有|A^{'}|=K\cdot|A|。

  1. 如果有两行完全相同，那么|A|=0。

运用性质1，我们有|A|=-|A^{'}|，由于这两行完全相同，我们又有A=A^{'}。

所以|A|=-|A|，所以|A|=0。

通过性质2推广，如果一行正好是另一行的K倍，那么有|A|=|A^{'}|

  1. 将其中一行加上另一行的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)。

应用

P6178 【模板】Matrix-Tree 定理

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;
}

P3317 [SDOI2014] 重建

题意简述：每条边概率连通，求最求概率最小生成树。

Solution

由上面的例题我们容易发现，矩阵树定理求的是对于图上所有的生成树T，设w(e)为e的边权，求解：

\sum_T\prod_{e\in T}w(e)

通过概率的加法/乘法原理，我们将所有的生成树的概率加起来，一颗生成树的概率为\prod_{e\in T}P_e\prod_{e\notin T}(1-P_e)，所以答案就是：

\sum_T\prod_{e\in T}P_e\prod_{e\notin T}(1-P_e)

推推式子，将式子化为\sum_T\prod_{e\in T}w(e)的样子。

\sum_T\prod_{e\in T}P_e\prod_{e\notin T}(1-P_e) =\sum_T\prod_{e\in T}P_e\frac{\prod_{e\in G}(1-P_e)}{\prod_{e\in T}(1-P_e)}

容易发现\prod_{e\in G}(1-P_e)是定值，将它移到外面。

=\prod_{e\in G}(1-P_e)\sum_T\frac{\prod_{e\in T}P_e}{\prod_{e\in T}(1-P_e)} =\prod_{e\in G}(1-P_e)\sum_T\prod_{e\in T}\frac{P_e}{(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;
}

P4336 [SHOI2016] 黑暗前的幻想乡

给定一张图，边有颜色，求正好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;
}