MatrixTree矩阵树定理

0. 前言

前置知识：行列式。

我已经 500 万年没写过这个笔记了，赶紧抽空写一个。

1. 生成树计数

1.1 无向图不带权

我们对于无向图，定义$D_{G}$ 表示图$G$ 的度数矩阵，即：

$$D_{G}(i,j)= \begin{cases} \text{deg}_{i} & (i=j) \\ \\ 0 & (i\neq j) \end{cases}$$

还有邻接矩阵，这里就不解释了，不过需要注意的是如果有重边的话应该算有多少条。

同时介绍 Kirchhoff 矩阵，指的对于一个图构造出来的一个矩阵。具体定义为度数矩阵减去邻接矩阵。

矩阵树定理说的是，一个图中的生成树个数等于其 Kirchhoff 矩阵的任意一个 代数余子式的行列式，说人话就是对其矩阵任意一个$n-1\times n-1$ 的子矩阵求行列式，注意如果直接对$n\times n$ 求行列式答案为$0$，因为一个图的 Kirchhoff 矩阵行列式为零，同时一个图的 Kirchhoff 矩阵的任一代数余子式的行列式相同。

模板题：SP104 HIGH - Highways - 洛谷，容斥题 P4336 [SHOI2016] 黑暗前的幻想乡 - 洛谷

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=620;
int n,m,a[MN][MN];

namespace HLS{

    int solve(){
        int ret=1,w=1;
        for(int i=1;i<=n;i++){
            for(int j=i+1;j<=n;j++){
                while(a[i][i]){
                    int div=a[j][i]/a[i][i];
                    for(int k=i;k<=n;k++){
                        a[j][k]=(a[j][k]-1ll*div*a[i][k]);
                    }
                    swap(a[i],a[j]);
                    w=-w;
                }
                swap(a[i],a[j]);
                w=-w;
            }
        }
        for(int i=1;i<=n;i++){
            ret=1ll*a[i][i]*ret;
        }
        ret=1ll*w*ret;
        return ret;
    }

}

void init(){
    for(int i=0;i<=n;i++){
        for(int j=0;j<=n;j++){
            a[i][j]=0;
        }
    }
}

void solve(){
    cin>>n>>m;
    init();
    for(int i=1;i<=m;i++){
        int u,v;
        cin>>u>>v;
        a[u][u]++;
        a[v][v]++;
        a[u][v]--;
        a[v][u]--;
    }
    n--;
    cout<<HLS::solve()<<'\n';
}

signed main(){
    int T;
    cin>>T;
    while(T--){
        solve();
    }
    return 0;
}

1.2 无向图带权

以下推广 Matrix-Tree 定理，推广到一般情况，Matrix Tree其实求的是所有可能生成树边权之积的和。由于每条生成树都算一次，权重是度数，分摊一下就是$1$，所以直接是生成树个数。

推广到带权图，假设图中边有权重$w_{e}$，我们希望求解加权生成树总和
，即对每棵生成树，把它包含的边权乘起来，再把所有生成树累加起来：

$$\sum\limits_{T}\prod_{e\in T} w_{e}$$

那么我们如何修改 Kirchhoff 矩阵呢，我们如下定义加权 Kirchhoff 矩阵：

• $(i,i)=\sum\limits_{(i,j)\in E} w(i,j)$。
• $(i,j)=-w(i,j)$，若无边默认为$0$。

然后同样的方法：去掉任意一行一列，求行列式，就可以得到结果。

P3317 [SDOI2014] 重建 - 洛谷

这道题是求刚好留下一棵生成树的概率，那么，即是求：

$$\begin{aligned} &\qquad \sum\limits_{T} \{\prod_{i\in T} p_{i}\prod_{i\notin T} (1-p_{i})\} \\ & = \sum\limits_{T} \{\prod_{i\in T} p_{i}\cdot \dfrac{\prod_{i}(1-p_{i})}{\prod_{i\in T}(1-p_{i})}\} \\ & =\prod_{i=1} (1-p_{i})\sum\limits_{T} \dfrac{p_{i}}{1-p_{i}} \end{aligned}$$

直接赋权即可，一个点的度数就是所有以它为端点的边的边权和。

#include<bits/stdc++.h>
#define db long double
using namespace std;
constexpr signed MN=55;
constexpr db eps=1e-8;
int n,m;
db a[MN][MN],ret=1;

namespace HLS{

    db solve(){
        db ret=1,w=1;
        for(signed i=1;i<=n;i++){
            for(signed j=i+1;j<=n;j++){
                    db div=a[j][i]/a[i][i];
                    for(signed k=i;k<=n;k++){
                        a[j][k]=(a[j][k]-div*a[i][k]);
                    }  
            }
        }
        for(signed i=1;i<=n;i++){
            ret=a[i][i]*ret;
        }
        return ret;
    }

}

signed main(){
    cin>>n;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            cin>>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;
            }
            if(fabs(1.0-a[i][j])<eps){
                a[i][j]=1.0-eps;
            }
            if(i<j) ret*=1.0-a[i][j];
            a[i][j]=a[i][j]/(1.0-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];
        }
    }
    n--;
    cout<<fixed<<setprecision(6)<<fabs(HLS::solve()*ret);
}

1.3 有向图版本

现在这个矩阵需要特殊定义，我们这里邻接矩阵的意义同有向图邻接矩阵，不妨设为$A$。

那么如果你要求的是外向树，那么矩阵就是：

$$M(i,j)= \begin{cases} \sum\limits_{k\neq i}w(k\to i) & (i=j) \\ \\ -w(j\to i) & (i\neq j) \end{cases}$$

即到该点的边权总和。

若是内向树，则为：

$$M(i,j)= \begin{cases} \sum\limits_{k\neq i}w(i\to k) & (i=j) \\ \\ -w(i\to j) & (i\neq j) \end{cases}$$

即从从该点出发的边权总和(出)。此外,既然是有向的,那么就需要指定根。在有向树的时候需要理解为指定根,结论是 : 去掉哪一行就是那一个元素为根。

例题：P4455 [CQOI2018] 社交网络 - 洛谷，有向图内向树 1 为根，代码如下：

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=620,MOD=1e4+7;
int n,m,t,a[MN][MN];

namespace HLS{

    int solve(){
        int ret=1,w=1;
        for(int i=1;i<=n;i++){
            for(int j=i+1;j<=n;j++){
                while(a[i][i]){
                    int div=a[j][i]/a[i][i];
                    for(int k=i;k<=n;k++){
                        a[j][k]=(a[j][k]-1ll*div*a[i][k]%MOD+MOD)%MOD;
                    }
                    swap(a[i],a[j]);
                    w=-w;
                }
                swap(a[i],a[j]);
                w=-w;
            }
        }
        for(int i=1;i<=n;i++){
            ret=1ll*a[i][i]*ret%MOD;
        }
        ret=1ll*w*ret;
        return (ret+MOD)%MOD;
    }

}

signed main(){
    cin>>n>>m;
    n--;
    for(int i=1;i<=m;i++){
        int u,v;
        cin>>u>>v;
        --u,--v;
        a[u][v]=(a[u][v]-1+MOD)%MOD;
        a[u][u]=(a[u][u]+1)%MOD;
    }
    cout<<HLS::solve();
    return 0;
}

1.4 无向带权与有向带权

P6178 【模板】Matrix-Tree 定理 - 洛谷

#include<bits/stdc++.h>
#define int long long
using namespace std;
constexpr int MN=620,MOD=1e9+7;
int n,m,t,a[MN][MN];

namespace HLS{

    int solve(){
        int ret=1,w=1;
        for(int i=1;i<=n;i++){
            for(int j=i+1;j<=n;j++){
                while(a[i][i]){
                    int div=a[j][i]/a[i][i];
                    for(int k=i;k<=n;k++){
                        a[j][k]=(a[j][k]-1ll*div*a[i][k]%MOD+MOD)%MOD;
                    }
                    swap(a[i],a[j]);
                    w=-w;
                }
                swap(a[i],a[j]);
                w=-w;
            }
        }
        for(int i=1;i<=n;i++){
            ret=1ll*a[i][i]*ret%MOD;
        }
        ret=1ll*w*ret;
        return (ret+MOD)%MOD;
    }

}

signed main(){
    cin>>n>>m>>t;
    n--;
    for(int i=1;i<=m;i++){
        int u,v,w;
        cin>>u>>v>>w;
        --u,--v;
        if(!t){
            a[u][v]=(a[u][v]-w+MOD)%MOD;
            a[v][u]=(a[v][u]-w+MOD)%MOD;
            a[u][u]=(a[u][u]+w+MOD)%MOD;
            a[v][v]=(a[v][v]+w+MOD)%MOD;
        }else{
            a[v][v]=(a[v][v]+w)%MOD;
            a[v][u]=(a[v][u]-w+MOD)%MOD;
        }
    }
    cout<<HLS::solve();
    return 0;
}