矩阵树定理学习笔记
用途
基尔霍夫矩阵树定理是在解决无向图的生成树计数问题的一种解决方案。
前置知识
主要就是行列式（及其性质）了。
行列式的定义
感性理解一下在这里就是一个矩阵换皮，并且包含了自己独特的求值方式。
这就是行列式写出来的形式了。可以看到，行列式就是把矩阵两边的大括号换成了一条竖杠罢了。
至于其特殊的求值方式，就是有一个求和式，对于一个(n imes n)的矩阵(M)，其行列式的值为：
行列式的性质
行列式有几个很有意思的性质，在后面实际应用的时候用的上。
-
互换一个矩阵的任意两行/列，其行列式值变为原值的相反数。
从上面的求值式子入手。我们很容易发现，交换了两行/列以后，两个值的绝对值不会变化，因为对于任意一个交换后的矩阵的行列式，所取到的(prodlimits_{i=1}^{n}M_{i,p_i})，都对应着原来两行交换的另一个(prodlimits_{i=1}^{n}M_{i,p_i})，这一对他们的(prod)里面包含的(M_{i,p_i})是同一个集合，只不过在顺序上有两个刚好换了而已。
但是，(prod)这一项不变，不代表前面的(-1^{ au(p)})（( au(p))表示(p)的逆序对数量）的不变。由上，我们可以想到把互换两行/列在证明中等价于一个排列交换掉两个数。而有一个很好证的结论，在(1)到(n)的全排列中，交换两个数所带来的逆序对数的变化是个奇数（可以通过考虑两个数之间与交换前后的两个数的大小关系来证，具体 留 作 读 者 自 证）。所以每次换行/列都会导致(-1)项的指数变化，即变号。
由此，我们还可以得到一个有趣的推论：
1.5. 矩阵有两行/列相同时，其行列式值为(0)：由于交换两行/列就会边号，但是交换这两行/列相同的行/列后行列式不变，所行列式的值只能为(0)。
-
向矩阵上的某一行/列上所有元素乘以一个数(k)，则行列式的值变为原值的(k)倍。
这个很好想，原式子里面的(prod)相等次数地包含了矩阵的所有行/列，所以可以从(prod)里面提出一个(k)就证毕了。
2.5. 矩阵中若有两行/列成比例，其行列式值为(0)：证明思路同推论1.5，只不过是先提出一行/列的关于另一行的倍数罢了。
-
（※重要）把矩阵上的一行/列加上另一行/列的(k)倍，其行列式值不变。
这个可以从求和式子的每一项的那一行的那个元素下手，把原来的求值式拆成两个式子：一个与原行列式相同，另一个所代表的矩阵中有两行成比例，值为0 。所以相比原来的行列式，值不变。
前置定义
- 度数矩阵(D)：一个(n imes n)的矩阵，其中(D_{i,i} = mathrm{degree}_i)，即(i)号点的度数，除此之外的矩阵元素为(0)。
- 邻接矩阵(A)：一个(n imes n)的矩阵，其中(A_{i,j} = ext{{直接链接i,j两点的边数/权值之类的}})。
- 基尔霍夫矩阵：一个(n imes n)的矩阵(K = D-A)。
基尔霍夫矩阵树定理
对于一个无向图(G)，其生成树数量等于其基尔霍夫矩阵减去任意一行及一列后，剩下的矩阵 所组成的行列式 的绝对值。
没错，就这么一句话。然而由于全网基本都没有给出过证明，唯一找到的一个证明在这，我也没看懂，所以用就完事了。
实际实现
我一开始兴高采烈的直接用定义去写然后发现全排列是(O(n!n))的，打完了才发现不行。
很容易发现定义的式子时间复杂度无法接受。所以我们要想一个办法快速求掉行列式的值。
首先由定义式，我们发现要想优化复杂度，肯定不能枚举全排列。然而全排列数量只随(n)而定。所以我们应该想些办法能让大部分的全排列后面的(prodlimits_{i=1}^{n}M_{i,p_i})值无用——即让他们为(0)。
考虑到刚才提到的行列式性质之一：
把矩阵上的一行/列加上另一行/列的(k)倍，其行列式值不变。
我们可以对矩阵做一下高斯消元，把它消成一个上三角矩阵。对于一个上三角矩阵，它所有的(prodlimits_{i=1}^{n}M_{i,p_i})中只有当(i = p_i)的时候才能取到非0的值（明显这个矩阵在(p_i<i)的元素都为(0)，然后如果有一个(i > p_i)，就必定有一个(i < p_i)，从而使积变为(0)）。所以说复杂度变为了高斯消元的(O(n^3))了。而且很明显，这个方法适用于所有的行列式，即它不是基尔霍夫矩阵的特殊情况。
此外，这种求生成树个数的题常常要取模，此时就需要在高斯消元的时候注意一下，浮点数是不可能的了。而且还要注意一下前面靠交换行而导致的变号。
例题1 洛谷P4336【SHOI2016】黑暗前的幻想乡
如果把每个公司的职权限制这一条件去掉，这就是一道求生成树个数的裸题，然而这道题加了这一条件，怎么办呢？
我们求的总个数也可以说是求了最多有(n-1)个公司参与修建的生成树，是肯定包含了只用了(n-2, n-3, dots)个公司的情况的。我们考虑到我们可以容易地求到最多有(n-1)家公司参与修建的方案数，是因为我们把所有公司的边都加进去了，那么只要我不加某些公司的边，然后直接矩阵树，不就是求了一定没有那些公司的方案数了嘛？这样子想，我们发现可以容易的求出最多有(x)家公司进行建设的方案数，最多->恰好这已经是人尽皆知的容斥套路了，最后我们求得的东西就是(sumlimits_{i=0}^{n-1}(-1)^{i+1}T(i))，其中(T(i))表示最多(i)个公司参加建设的方案数。
此外，这道题由于公司的职权可能有重复，所以有重边。这说明了矩阵树定理可以有重边（而且很多时候所谓树边权值也能化为重边考虑）（其实带重边/权值的矩阵树定理就是变元矩阵树定理了）。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
namespace ztd{
using namespace std;
const int mod = 1e9+7;
typedef long long ll;
template<typename T> inline T read(T& t) {//fast read
t=0;short f=1;char ch=getchar();double d = 0.1;
while (ch<'0'||ch>'9') {if (ch=='-') f=-f;ch=getchar();}
while (ch>='0'&&ch<='9') t=t*10+ch-'0',ch=getchar();
if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
t*=f;
return t;
}
inline ll ksm(ll x, ll y){
ll ret = 1;
while(y){if(y & 1) ret = ret * x % mod; y >>= 1; x = x*x%mod;}
return ret;
}
}
using namespace ztd;
const int maxn = 25;
int n, m[maxn], ip[maxn][1007][2];
struct matrix{//以前自己封装的矩阵板子
long long a[25][25], n, m;
matrix(){ memset(a,0,sizeof(a));} //初始化
inline void init(int x, int y){memset(a,0,sizeof(a)); n = x, m = y;}
matrix operator - (const matrix &Y)const{
matrix X = *this;
for(int i = 1; i <= X.n; ++i){
for(int j = 1; j <= X.m; ++j){
X.a[i][j] = (X.a[i][j] - Y.a[i][j] + mod) % mod;
}
}
return X;
}
inline int Gauss(){//高斯消元
int sign = 1;//记录变号的次数
for(int i = 1; i < n; ++i){//要去掉一行，所以只到n-1
if(!a[i][i]) {//由于要构造上三角，所以尝试找到一行和现在这一行交换，
//由前面的性质，行列式值会变成原相反数。
for(int j = i+1; j < n; ++j){
if(!a[j][i]) continue;
sign = -sign;
swap(a[i], a[j]);
break;
}
}
if(!a[i][i]) break;//如果有一个a[i][i]在更新完之后还是0，那么积就为0.
ll inv = ksm(a[i][i], mod-2);
for(int j = i+1; j < n; ++j){//消掉底下的a[i][i]
ll tmp = a[j][i] * inv % mod;
for(int k = i; k < n; ++k)
a[j][k] = (a[j][k] - (a[i][k]*tmp%mod) + mod) % mod;
}
}
return sign;
}
inline int get_det(){//获求行列式的值
int sign = Gauss();//获取一下变号情况
ll ret = 1;
for(int i = 1; i < n; ++i)
ret = ret * a[i][i] % mod;
return ((ret*sign)+mod)%mod;
}
}deg, edg, Kirchhoff;
ll ans;
void dfs(int x, int now){ //dfs每一个公司的参加情况
if(x == n) return;
//不选
dfs(x+1, now);
//选
for(int i = 1, xx, yy; i <= m[x]; ++i){
xx = ip[x][i][0], yy = ip[x][i][1];
++deg.a[xx][xx]; ++deg.a[yy][yy];
++edg.a[xx][yy]; ++edg.a[yy][xx];
}
Kirchhoff = deg-edg;//K = D-A， 整个矩阵树定理的核心，就是一行代码。
ans = (ans + ((n-now)&1 ? -1:1) * Kirchhoff.get_det() + mod) % mod;
dfs(x+1, now+1);
//回溯
for(int i = 1, xx, yy; i <= m[x]; ++i){
xx = ip[x][i][0], yy = ip[x][i][1];
--deg.a[xx][xx]; --deg.a[yy][yy];
--edg.a[xx][yy]; --edg.a[yy][xx];
}
}
signed main(){
read(n); deg.init(n,n); edg.init(n,n);
for(int i = 1; i < n; ++i) {
read(m[i]);
for(int j = 1; j <= m[i]; ++j)
read(ip[i][j][0]); read(ip[i][j][1]);
}
dfs(1, 0);
cout << ans << '
';
return 0;
}
例题2 洛谷P6624【省选联考 2020 A 卷】作业题
这道题从他给的那个要求的式子来看，就十分的……缝合。
右边是个GCD。可以用欧拉反演变成(sumlimits_{d|gcd}varphi(d))，可以线性筛后(O(1))求值。这不是我们的重点。
把视角转向左边半边的式子。我们想起了我们前面求用矩阵树定理解的那些东西（如果把重边视做权值，毕竟这题没有重边）是(sumlimits_{T}prodlimits_{i=1}^{n-1}w_{e_i})，而这道题求的是(sumlimits_{T}sumlimits_{i=1}^{n-1}w_{e_i})。有什么方法可以解决这个问题呢？
有！这里有一个技巧，就是把所有边边权变为一个一次生成函数(w_ix+1)，最后取答案的一次项。这样子两个多项式相乘，他们的一次项就是相加，达到了我们想要的效果了。而且由于只取答案一次项，所以维护的”多项式“只用维护一次项和常数项即可。这些一次“多项式”的加减乘比较好些，除法可以用多项式求逆搞一下，就能推出式子了。其他的就是封装好，矩阵树定理完事。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
const int maxn = 33, maxn2 = 1007, maxm = 200007;
namespace ztd{
using namespace std;
typedef long long ll;
template<typename T> inline T read(T& t) {//fast read
t=0;short f=1;char ch=getchar();double d = 0.1;
while (ch<'0'||ch>'9') {if (ch=='-') f=-f;ch=getchar();}
while (ch>='0'&&ch<='9') t=t*10+ch-'0',ch=getchar();
if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
t*=f;
return t;
}
const int mod = 998244353;
inline ll ksm(ll x, ll y){
ll ret = 1;
while(y){if(y & 1) ret = ret * x % mod; y >>= 1; x = x*x%mod;}
return ret;
}
}
using namespace ztd;
namespace NumberTheroy{//预处理phi
int phi[maxm], fct[maxm], prime[maxm], pcnt;
bool vis[maxm];
inline void add_factor(int x){//统计各个数的因子
for(int i = 1; i*i <= x; ++i){
if(x%i == 0){
++fct[i];
if(i*i != x) ++fct[x/i];
}
}
}
inline void preprocess(int uplimit){//线性筛
vis[1] = phi[1] = 1;
for(int i = 2; i <= uplimit; ++i){
if(!vis[i]) phi[i] = i-1, prime[++pcnt] = i;
for(int j = 1; j <= pcnt && prime[j]*i <= uplimit; ++j){
vis[i*prime[j]] = 1;
if(i % prime[j] == 0){
phi[prime[j]*i] = phi[i]*prime[j];
break;
}
phi[prime[j]*i] = phi[i]*(prime[j]-1);
}
}
}
}
using namespace NumberTheroy;
int n, m, x[maxn2], y[maxn2], w[maxn2];
struct node{//一次多项式
ll x, y;
//常数项⬆ ⬆一次项
node operator +(const node a){
return (node){(x+a.x)%mod, (y+a.y)%mod};
}
node operator -(const node a){
return (node){(x-a.x+mod)%mod, (y-a.y+mod)%mod};
}
node operator *(const node a){
return (node){x*a.x%mod , ((x*a.y%mod) + (y*a.x%mod) )%mod};
}
node operator *(const int a){
return (node){(a*x+mod)%mod, (a*y+mod)%mod};
}
node operator /(const node a){//多项式求逆
int inv = ksm(a.x, mod-2);
return (node){x*inv%mod, (y*a.x%mod - x*a.y%mod + mod)%mod *inv%mod *inv%mod};
}
};
struct matrix{
node a[maxn][maxn]; int n, m;
matrix(){ memset(a,0,sizeof(a));} //初始化
inline void init(int x, int y){memset(a,0,sizeof(a)); n = x, m = y;}
inline int Gauss(){//高斯消元，和上面没什么差别，就是原来的int的矩阵变成了多项式矩阵。
int sign = 1;
for(int i = 1; i < n; ++i){
if(!a[i][i].x){
for(int j = i+1; j < n; ++j){
if(!a[j][i].x) continue;
sign = -sign;
for(int k = 1; k < n; ++k) swap(a[i][k], a[j][k]);
break;
}
}
if(!a[i][i].x) break;
node inv = ((node){1,0}) / a[i][i];
for(int j = i+1; j < n; ++j){
node tmp = a[j][i] * inv;
for(int k = i; k < n; ++k) a[j][k] = a[j][k] - a[i][k]*tmp;
}
}
return sign;
}
inline node get_det(){
int sign = Gauss(); node ret = (node){1,0};
for(int i = 1; i < n; ++i) ret = ret * a[i][i];
return ret*sign;
}
}G;
ll ans;
inline ll solve(int X){
G.init(n,n);
for(int i = 1; i <= m; ++i){
if(w[i] % X) continue;
G.a[x[i]][x[i]] = G.a[x[i]][x[i]] + (node){1, w[i]};
G.a[y[i]][y[i]] = G.a[y[i]][y[i]] + (node){1, w[i]};
G.a[x[i]][y[i]] = G.a[x[i]][y[i]] - (node){1, w[i]};
G.a[y[i]][x[i]] = G.a[y[i]][x[i]] - (node){1, w[i]};
}
node ret = G.get_det();
return ret.y;
}
signed main(){
read(n); read(m);
int mxw = 0; ll ans = 0;
for(int i = 1; i <= m; ++i){
read(x[i]); read(y[i]); read(w[i]);
add_factor(w[i]);
mxw = max(mxw, w[i]);
}
preprocess(mxw);
for(int i = 1; i <= mxw; ++i) if(fct[i] >= n-1){
ans = (ans + 1ll * phi[i] * solve(i) % mod) % mod;
}
cout << ans << '
';
return 0;
}
矩阵树定理在有向图的扩展
前面讲的都是无向图，其实我们还可以在有向图中用矩阵树定理求出生成“树”的数量。
有向图的生成树
首先明确一点，有向图没有生成树。所以我们求得就是它的生成外向树/生成内向树。其中：
- 内向树表示把树上的所有边替换成指向父亲的有向边所形成的树
- 外向树表示把树上的所有边替换成指向所有儿子的有向边所形成的树
很明显，这一类“树”是要明确根为哪一个点的。
前置定义
- 入度/出度矩阵(D_{in} / D_{out})：同原定理的度数矩阵(D)，只不过对于每一个(D_{i,i})记录的变为第(i)个点的入度/出度了。
- 基尔霍夫矩阵的扩展：(K_{in} = D_{in} - A, K_{out} = D_{out} - A) 。
定理内容
对于一个有向图(G)，
- 其根为(rt)的生成外向树数量等于其(K_{in})减去第(rt)行及第(rt)列后，剩下的矩阵 所组成的行列式 的绝对值;
- 其根为(rt)的生成内向树数量等于其(K_{out})减去第(rt)行及第(rt)列后，剩下的矩阵 所组成的行列式 的绝对值。
例题3 P6178 【模板】Matrix-Tree 定理
一个有边权、自环的有向图，求生成外向树。那就是上面的模板了。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
namespace ztd{
using namespace std;
typedef long long ll;
const int mod = 1e9+7;
template<typename T> inline T read(T& t) {//fast read
t=0;short f=1;char ch=getchar();double d = 0.1;
while (ch<'0'||ch>'9') {if (ch=='-') f=-f;ch=getchar();}
while (ch>='0'&&ch<='9') t=t*10+ch-'0',ch=getchar();
if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
t*=f;
return t;
}
inline ll ksm(ll x, ll y){
ll ret = 1;
while(y){if(y & 1) ret = ret * x % mod; y >>= 1; x = x*x%mod;}
return ret;
}
}
using namespace ztd;
const int maxn = 307;
int n, m, t;
struct matrix{
ll a[maxn][maxn], n;
inline void init(int nn){memset(a,0,sizeof(a)); n = nn;}
inline int Gauss(){
int sign = 1;
for(int i = 2; i <= n; ++i){
if(!a[i][i]) {
for(int j = i+1; j <= n; ++j){
if(!a[j][i]) continue;
sign = -sign;
swap(a[i], a[j]);
break;
}
}
if(!a[i][i]) break;
ll inv = ksm(a[i][i], mod-2);
for(int j = i+1; j <= n; ++j){
ll tmp = a[j][i] * inv % mod;
for(int k = i; k <= n; ++k)
a[j][k] = (a[j][k] - (a[i][k]*tmp%mod) + mod) % mod;
}
}
return sign;
}
inline ll get_det(){
int sign = Gauss();
ll ret = 1;
for(int i = 2; i <= n; ++i) ret = ret*a[i][i]%mod;
return ((ret*sign)+mod) % mod;
}
}K;
signed main(){
read(n); read(m); read(t); K.init(n);
int xx, yy, ww;
for(int i = 1; i <= m; ++i){
read(xx); read(yy); read(ww);
K.a[yy][yy] = (K.a[yy][yy] + ww) % mod; K.a[xx][yy] = (K.a[xx][yy]-ww+mod)%mod;
if(!t) K.a[xx][xx] = (K.a[xx][xx] + ww) % mod, K.a[yy][xx] = (K.a[yy][xx]-ww+mod)%mod;
}
cout << (K.get_det()) << '
';
return 0;
}
其他题目
- P4455 【CQOI2018】社交网络
一道有向扩展的板子题，相比上面的差异就在于求的是生成内向树 - 未完待续