Matrix – Tree 定理学习笔记

ooliver 发布于 22 小时前 130 次阅读 OI

AI 摘要

一个行列式，竟能数清一张图的所有生成树？矩阵树定理，数学的魔法开关，带权、有向、内向、外向——统统化为行列式。从生成树计数到MST统计，看完这篇笔记，你也将掌握这门“降维打击”的技巧。

前置芝士

Matrix - Tree 定理，中文名：矩阵树定理。

度数矩阵

对于无向图，度数矩阵 $D$ 的定义是：

$$
D_{i,j}=
\begin{cases}
deg_i &\text{if } i=j \\
0 &\text{if } i\not =j \\
\end{cases}
$$

对于有向图，分入度矩阵 $D^{in}$ 和出度矩阵 $D^{out}$，定义如下：

$$
D^{in}_{i,j}=
\begin{cases}
in_i &\text{if } i=j \\
0 &\text{if } i\not =j \\
\end{cases}
$$

$$
D^{out}_{i,j}=
\begin{cases}
out_i &\text{if } i=j \\
0 &\text{if } i\not =j \\
\end{cases}
$$

邻接矩阵

对于无向图，邻接矩阵 $A$ 的定义为：

$$
A_{i,j}=A_{j,i}=\# e(i,j)
$$

其中 $\# e(i,j)$ 表示 $i,j$ 两点之间所连边数。

对于有向图，邻接矩阵 $A$ 的定义为：

$$
A_{i,j}=\# e(i,j)
$$

其中 $\# e(i,j)$ 表示从点 $i$ 指向 $j$ 的边数。

拉普拉斯矩阵

说到拉普拉斯就不得不提到一位恶魔系拉普拉斯妖学妹

对于无向图，拉普拉斯矩阵 $L$ 的定义为：

$$
L=D-A
$$

对于有向图，入度拉普拉斯矩阵 $L^{in}$ 和出度拉普拉斯矩阵 $L^{out}$ 的定义为：

$$
L^{in}=D^{in}-A
$$

$$
L^{out}=D^{out}-A
$$

Matrix - Tree 定理

对于无向图，其生成树的个数为（$k\in [1,n]$)：

$$
\text{det }L_{[n] \setminus \{k\},[n] \setminus \{k\}}
$$

其中 $\text{det }$ 表示行列式求值，$L_{[n] \setminus \{k\},[n] \setminus \{k\}}$ 表示将矩阵 $L$ 删除第 $k$ 行与第 $k$ 列后得到的矩阵。

对于有向图，其以 $k$ 为根的内向生成树（根向生成树）的个数为：

$$
\text{det }L^{out}_{[n] \setminus \{k\},[n] \setminus \{k\}}
$$

其以 $k$ 为根的外向生成树（叶向生成树）的个数为：

$$
\text{det }L^{in}_{[n] \setminus \{k\},[n] \setminus \{k\}}
$$

带权拓展

问题：对于一个带权图，定义一个生成树的权值为其所有边权之积，求其所有生成树权值之和。

这里我们修改一下度数矩阵的定义即可，不在存储边数，而是所有边边权之和，按照这种方式得到拉普拉斯矩阵后其余和上面的式子是一模一样的。

例题

P6178 【模板】Matrix-Tree 定理

带权拓展的模板。注意，这里有向图生成树的定义是以 $1$ 为根的外向树，所以我们要看入度，并删除拉普拉斯矩阵的第一行、第一列。

代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define int long long
const int N=305;
const int mod=1e9+7;
int gh[N][N];

int det(int n){
    int res=1,f=1;
    for(int i=1;i<=n;i++){
        for(int j=i+1;j<=n;j++){
            while(gh[j][i]){
                int t=gh[i][i]/gh[j][i];
                for(int k=i;k<=n;k++) (((gh[i][k]-=gh[j][k]*t)%=mod)+=mod)%=mod;
                swap(gh[i],gh[j]);
                f=-f;
            }
        }
        if(!gh[i][i]) return 0;
        (res*=gh[i][i])%=mod;
    }
    return ((res*f)%mod+mod)%mod;
}

signed main(){
    int n,m,t;
    cin>>n>>m>>t;
    for(int i=1,u,v,w;i<=m;i++){
        cin>>u>>v>>w;
        u--,v--;
        if(t) (gh[v][v]+=w)%=mod,(((gh[u][v]-=w)%=mod)+=mod)%=mod;
        else (gh[v][v]+=w)%=mod,(((gh[u][v]-=w)%=mod)+=mod)%=mod,
             (gh[u][u]+=w)%=mod,(((gh[v][u]-=w)%=mod)+=mod)%=mod;
    }
    cout<<det(n-1);
    return 0;
}

P3317 [SDOI2014] 重建

推一下式子：

$$
\begin{aligned}
P&=\sum_{T \text{ 是生成树}} \left( \prod_{e \in T} p_e \prod_{e \notin T}(1 - p_e)\right) \\
&=\left( \prod_{e \in E}(1 - p_e) \right) \cdot \sum_{T \text{ 是生成树}} \left( \prod_{e \in T}\frac{p_e}{1 - p_e}\right)
\end{aligned}
$$

后面就是矩阵树定理的板子了，代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define double long double
const int N=105;
const double eps=1e-9;
double gh[N][N],res=1,f=1;

void det(int n){
    for(int i=1;i<=n;i++){
        int now=i;
        for(int j=i+1;j<=n;j++) if(fabs(gh[j][i])>fabs(gh[now][i])) now=j;
        swap(gh[i],gh[now]);
        if(i!=now) f=-f;
        if(gh[i][i]<eps){
            res=0;
            return;
        }
        for(int j=i+1;j<=n;j++){
            if(gh[j][i]==0) continue;
            double t=gh[j][i]/gh[i][i];
            for(int k=i;k<=n;k++) gh[j][k]-=t*gh[i][k];
        }
        res*=gh[i][i];
    }
}

signed main(){
    int n;
    cin>>n;
    for(int u=1;u<=n;u++) for(int v=1;v<=n;v++){
        double w;
        cin>>w;
        if(w>1-eps) w=1-eps;
        gh[v][v]+=w/(1-w),gh[u][v]-=w/(1-w);
        if(u>v) res*=(1-w);
    }
    det(n-1);
    cout<<fixed<<setprecision(5)<<res*f;
    return 0;
}

P4455 [CQOI2018] 社交网络

外向生成树计数板子题。

代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define int long long
const int N=305;
const int mod=1e4+7;
int gh[N][N];

int det(int n){
    int res=1,f=1;
    for(int i=1;i<=n;i++){
        for(int j=i+1;j<=n;j++){
            while(gh[j][i]){
                int t=gh[i][i]/gh[j][i];
                for(int k=i;k<=n;k++) (((gh[i][k]-=gh[j][k]*t)%=mod)+=mod)%=mod;
                swap(gh[i],gh[j]);
                f=-f;
            }
        }
        if(!gh[i][i]) return 0;
        (res*=gh[i][i])%=mod;
    }
    return ((res*f)%mod+mod)%mod;
}

signed main(){
    int n,m;
    cin>>n>>m;
    for(int i=1,u,v,w;i<=m;i++){
        cin>>v>>u;
        u--,v--;
        (gh[v][v]+=1)%=mod,(((gh[u][v]-=1)%=mod)+=mod)%=mod;
    }
    cout<<det(n-1);
    return 0;
}

P4336 [SHOI2016] 黑暗前的幻想乡

考虑容斥原理：

$$
所有公司都没有被漏掉的生成树方案数=\sum_{S\subseteq U} (-1)^{|S|} \times (S\ 内的公司被漏掉的生成树方案数)
$$

之后就是矩阵树定理了。

代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define int long long
const int N=20;
const int mod=1e9+7;
vector<pair<int,int>> ve[N];
int gh[N][N],ans;

int det(int n){
    int res=1,f=1;
    for(int i=1;i<=n;i++){
        for(int j=i+1;j<=n;j++){
            while(gh[j][i]){
                int t=gh[i][i]/gh[j][i];
                for(int k=i;k<=n;k++) (((gh[i][k]-=gh[j][k]*t)%=mod)+=mod)%=mod;
                swap(gh[i],gh[j]);
                f=-f;
            }
        }
        if(!gh[i][i]) return 0;
        (res*=gh[i][i])%=mod;
    }
    return ((res*f)%mod+mod)%mod;
}

signed main(){
    int n;
    cin>>n;
    for(int i=1;i<n;i++){
        int m;
        cin>>m;
        for(int j=1;j<=m;j++){
            int u,v;
            cin>>u>>v;
            ve[i].push_back({u,v});
        }
    }
    for(int s=0;s<(1<<(n-1));s++){
        memset(gh,0,sizeof gh);
        int cnt=0;
        for(int i=0;i<n-1;i++){
            if(!((s>>i)&1)) continue;
            cnt++;
            for(auto[u,v]:ve[i+1]) (gh[u][u]+=1)%=mod,(gh[v][v]+=1)%=mod,(gh[u][v]+=mod-1)%=mod,(gh[v][u]+=mod-1)%=mod;
        }
        ans=((ans+((n-cnt-1)&1?-1:1)*det(n-1))%mod+mod)%mod;
    }
    cout<<ans;
    return 0;
}

P2143 [JSOI2010] 巨额奖金 & P4208 [JSOI2008] 最小生成树计数

考虑所有的最小生成树的边权序列排序后都是一致的，所以我们先跑一遍最小生成树，对于这个最小生成树的每一种边权 $w$，我们在新图中先把原最小生成树中边权不为 $w$ 的边上的点进行缩点，再按照原图中所有边权为 $w$ 的边建图，跑矩阵树求生成树个数即可。

代码：

C++
#include<bits/stdc++.h>
using namespace std;

#define int long long
const int N=105,M=1005;
const int mod=31011;
int n,m,ans=1;
int gh[N][N],nd[N];
map<int,bool> mp;

struct edge{
    int u,v,w;
    bool operator<(const edge b)const{
        return w<b.w;
    }
}e[M];

vector<edge> mst;

struct DSU{
    int fa[N],siz[N];
    void init(){
        for(int i=1;i<=n;i++) fa[i]=i,siz[i]=1;
    }
    int find(int x){
        while(x!=fa[x]) x=fa[x]=fa[fa[x]];
        return x;
    }
    bool merge(int x,int y){
        int fx=find(x),fy=find(y);
        if(fx==fy) return 0;
        if(siz[fx]>siz[fy]) swap(fx,fy);
        return fa[fy]=fx,1;
    }
}dsu;

int det(int n){
    int res=1,f=1;
    for(int i=1;i<=n;i++){
        for(int j=i+1;j<=n;j++){
            while(gh[j][i]){
                int t=gh[i][i]/gh[j][i];
                for(int k=i;k<=n;k++) (((gh[i][k]-=gh[j][k]*t)%=mod)+=mod)%=mod;
                swap(gh[i],gh[j]);
                f=-f;
            }
        }
        if(!gh[i][i]) return 0;
        (res*=gh[i][i])%=mod;
    }
    return ((res*f)%mod+mod)%mod;
}

bool kk(){
    int res=0;
    dsu.init();
    sort(e+1,e+1+m);
    for(int i=1;i<=m;i++){
        if(!dsu.merge(e[i].u,e[i].v)) continue;
        mp[e[i].w]=1,res++;
        mst.push_back(e[i]);
        if(res==n-1) return 1;
    }
    return 0;
}

signed main(){
    cin>>n>>m;
    for(int i=1;i<=m;i++) cin>>e[i].u>>e[i].v>>e[i].w;
    if(!kk()){
        cout<<0;
        return 0;
    }
    for(auto[val,tag]:mp){
        int idx=0;
        dsu.init();
        memset(gh,0,sizeof gh);
        for(auto[u,v,w]:mst){
            if(w==val) continue;
            dsu.merge(u,v);
        }
        for(int i=1;i<=n;i++) if(dsu.find(i)==i) nd[i]=++idx;
        for(int i=1;i<=n;i++) if(dsu.find(i)!=i) nd[i]=nd[dsu.find(i)];
        for(int i=1;i<=m;i++){
            if(e[i].w!=val) continue;
            int x=nd[e[i].u],y=nd[e[i].v];
            (gh[x][x]+=1)%=mod,(gh[y][y]+=1)%=mod,
            (gh[x][y]+=mod-1)%=mod,(gh[y][x]+=mod-1)%=mod;
        }
        (((ans*=det(idx-1))%=mod)+=mod)%=mod;
    }
    cout<<ans;
    return 0;
}