@算法 - 7@ matrix - tree 定理（矩阵树定理）

标签：算法   fabs   put   end   push   有一个   了解   %s   答案   

目录

• @0 - 参考资料@
• @0.5 - 你所需要了解的线性代数知识@
• @1 - 定理主体@
  • @证明 part - 1@
  • @证明 part - 2@
  • @证明 part - 3@
  • ＠证明 part - 4＠
• @2 - 一些简单的推广@
• ＠3 - 例题与应用＠

@0 - 参考资料@

MoebiusMeow 的讲解（超喜欢这个博主的！）
网上找的另外一篇讲解

@0.5 - 你所需要了解的线性代数知识@

什么是矩阵？
什么是高斯消元？这个虽然与主题无关，但是求解行列式需要用。
什么是行列式？行列式采用排列的方式定义。但是为了简洁（我懒得写），我们直接使用拉普拉斯展开。

设一个 n*n 矩阵 A：
\[A = \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}\]
对第 i 行（i 可以是任意的）进行展开，得到 A 的行列式 det(A)为：
\[\det(A) = \sum_{j=1}^{n}(-1)^{i+j}a_{ij}M_{ij}\]
其中 \(M_{ij}\) 表示去掉第 i 行第 j 列后剩下的矩阵的行列式，称为 \(a_{ij}\) 的余子式。

行列式的一些性质（我真的懒得写证明）：
（1）矩阵 A 的行列式等于其转置的行列式（可以理解为行和列是等价的）。
（2）把矩阵 A 的某一行变为 k 倍，行列式变为 k 倍。
（3）把矩阵 A 的第 i 行加第 j 行的若干倍，行列式不变。
（4）交换两行，行列式变号（正变负，负变正）。
（5）上/下三角矩阵的行列式等于主对角线元素乘积。

（1）表示行列式中行和列是等价的。
（2）、（3）和（4）都是关于初等行变换的。
（5）提供了一种计算特殊矩阵行列式的方法。

因此，我们可以通过高斯消元将一个矩阵消成三角矩阵 O(n^3) 求行列式。

@1 - 定理主体@

先考虑最简单的：无重边无自环的无向图 G。（简称三无图？）

我们定义这个无向图 G 的度数矩阵 \(D\)。只有 \(D_{ii}\) 为点 i 的度数，其他项为 0。

我们定义这个无向图 G 的邻接矩阵 \(A\)。\(A_{ij}\) 为 1 当且仅当点 i 与点 j 有连边。

我们定义这个无向图 G 的基尔霍夫矩阵 \(L\)。\(L = D - A\)。

matrix - tree 定理：这个无向图的生成树个数等于任意一个主对角线上的元素的余子式。

@证明 part - 1@

我们先来证明一个简单的引理：任意一个基尔霍夫矩阵的行列式总等于零。

观察到基尔霍夫矩阵，其每一行每一列的和等于零。

因此我们可以将其他行加到第一行，第一行就变成全零的。

再把第一行展开就可以证明了。

@证明 part - 2@

再证明这样一个引理：如果 G 不连通，主对角线上的余子式总等于零。

我们将基尔霍夫矩阵 L 作这样的变换：将每一个连通块的点都放在一起。
每次交换 i, j 两行，再交换 i, j 两列，就可以将点 i 与点 j 的位置对换，且行列式两次变号后还是保持不变。

这样最终基尔霍夫矩阵长这样：
\[L = \begin{bmatrix}A_1 & 0 & \cdots & 0\\ 0 & A_2 & \cdots & 0 \\ \vdots & \vdots& \ddots & \vdots \\ 0 & 0 & \cdots & A_k \end{bmatrix}\]

这个时候，\(\det(L)=\det(A_1)\det(A_2)\dots\det(A_k)\)。至于为什么直观感受一下，对角线的乘积嘛。
其实证明也比较简单，将所有 A 高斯消元消成三角矩阵就可以了。

因为不连通，所以 \(2 \leq k\)。
我们的余子式中至少包含一个完整的 A 矩阵。
而每一个矩阵 A 都是基尔霍夫矩阵。所以最终的余子式一定等于零。

@证明 part - 3@

接下来我们证明最后一个引理：如果 G 是一棵树，主对角线上的余子式总等于一。

我们先把要展开的那一行的点作为根结点，然后以根为起点进行 bfs。
按照 bfs 遍历的逆序将结点排列（最后遍历的放最前面）。
一样还是每次交换 i, j 两行再交换 i, j 两列交换 i, j 两点的位置。

然后，我们从左往右扫描每一列。每次扫描到某一列，我们将这一列加到这个结点的父亲那一列上（如果它的父亲是根结点就什么也不做）。

通过归纳法可以证明，扫描到第 i 列的时候，第 i 列的第 1~i-1 行全是 0，第 i 行为 1，第 i+1~n 行上除了它父亲那一行为 -1 其他都为 0。
直观理解也可以理解。

于是最后剩下的矩阵是三角矩阵且主对角线上都是 1。
行列式自然就等于 1 了。

＠证明 part - 4＠

要证明矩阵树定理，我们要先知道 Cauchy - Binet定理 ：
对于 n*m 矩阵 A 与 m*n 矩阵 B （m < n）：
\[\det(AB)=\sum_{|s|=m}\det(A_{s*})\det(B_{s*})=\sum_{|s|=m}\det(A_{s*}*B_{s*})\]
表示我们选择 A 中的 m 行构成的 m*m 的方阵乘上选择 B 中相应的 m 列构成的 m*m 的方阵之和等于 AB 的行列式。

该定理的证明超出了我的理解范围。

然后我们开始我们正式的证明：
我们构造 m*n 矩阵 \(B\)，其中 m 是边数，n 是点数。
如果第 i 条边为 \((u_i, v_i)\)，矩阵中就有 \(b_{iu_i}=1, b_{iv_i}=-1\)（因为这条边是无向边，所以哪一个是 -1 哪一个是 1 没有什么实际影响，只要保证一个为 1 另一个为 -1 就可以了。）

然后，有 \(L = B^TB\)。为什么呢？可以把矩阵乘法的式子列出来看。

再根据上面那个定理，就有：
\[M_{ii}=\det(B^T_iB_i)=\sum_{|s|=n-1}\det((B_i^T)_{s*}(B_i)_{*s})\]

右边那个其实就是从 m 条边里面选择 n-1 条边构成的图，展开的余子式。根据我们的几个引理，只有当这个图是树的时候，会对答案产生贡献且贡献为 1。所以余子式等于生成树个数。

@2 - 一些简单的推广@

首先最关注的应该是重边自环。
根据上面的推导过程，重边是没有影响。
自环的话，根据我们构造的 B 矩阵，你把 Bij++ 又把 Bij-- 过后，实际上这个自环就不会造成任何影响了。但是实现上可能要特殊处理一下。

然后就是有向图的推广。有向图的生成树，我们这里只考虑以结点 i 为根的外向生成树。其他情况可以类比得出。
考虑我们证明树的余子式总为 1 时，我们是用一个点的度数去消去它父亲的邻接矩阵。因此，我们将邻接矩阵改成只记出边的矩阵，将度数矩阵改为只记入度的矩阵（当然因为矩阵的行列是等价的所以你也不需要特别去弄清楚谁是出边谁是入边，多试试就可以了）。
注意我们如果是有向图，不能随便展开行，展开的那一行就是根。

再然后，就是其他类型的生成树统计而不仅仅是个数统计。最基础的应该是每条边带边权，求不同情况生成树的边权乘积的总和。可以发现边权为 1 时就是我们的正常个数统计，因此我们可以直接将该加 1 减 1 的地方换成边权即可。

还有一个高斯消元上的细节问题。因为我们要求解的是行列式，所以用不着使用将某一行变 k 倍的初等行变换，让每一行的第一个非零元变成 1。
因此就有一个问题，是否我们可以在不一定存在逆元的模数（即非质数）意义下作行列式的求解。
我们知道使用辗转相除法，在 log 次就可以使得某一个数为 0。
因此我们对两行使用辗转相除法，在 O(n^3log n) 的时间复杂度实现非质数模数的行列式求解。

＠3 - 例题与应用＠

一开始竟然把 gauss 拼写成了 guess，果然太久没有写过这东西了……。

首先是无向图的生成树计数模板题（可能建图不是那么模板？）。
bzoj - 4031 : [HEOI2015]小Z的房间

#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN = 9;
const int MAXV = 81;
const int MOD = int(1E9);
const int dx[] = {0, 0, 1, -1};
const int dy[] = {1, -1, 0, 0};
int mat[MAXV][MAXV];
int gauss(int r, int c) {
    int nwr = 0, nwc = 0, f = 1;
    for(int i=0;i<r;i++)
        for(int j=0;j<c;j++)
            mat[i][j] = (mat[i][j] + MOD)%MOD;
    while( nwr < r && nwc < c ) {
/*
        for(int i=0;i<r;i++)
            for(int j=0;j<c;j++)
                printf("%9d%c", mat[i][j], (j + 1 == c) ? ‘\n‘ : ‘ ‘);
        puts("");
*/
        int mxr = nwr;
        for(int i=nwr+1;i<r;i++)
            if( mat[i][nwc] > mat[mxr][nwc] ) mxr = i;
        if( mat[mxr][nwc] ) {
            if( mxr != nwr ) {
                for(int i=nwc;i<c;i++)
                    swap(mat[nwr][i], mat[mxr][i]);
                f *= -1;
            }
            for(int i=nwr+1;i<r;i++)
                while( mat[i][nwc] ) {
                    if( mat[nwr][nwc] < mat[i][nwc] ) {
                        for(int j=nwc;j<c;j++)
                            swap(mat[nwr][j], mat[i][j]);
                        f *= -1; continue;
                    }
                    else {
                        int x = mat[nwr][nwc]/mat[i][nwc];
                        for(int j=nwc;j<c;j++)
                            mat[nwr][j] = (mat[nwr][j] + MOD - 1LL*x*mat[i][j]%MOD)%MOD;
                    }
                }
            nwr++;
        }
        nwc++;
    }
    int ans = 1;
    for(int i=0;i<r;i++)
        ans = 1LL*ans*mat[i][i]%MOD;
    return (ans*f + MOD)%MOD;
}
int num[MAXN][MAXN];
char str[MAXN];
int main() {
    int n, m, cnt = 0;
    scanf("%d%d", &n, &m);
    for(int i=0;i<n;i++) {
        scanf("%s", str);
        for(int j=0;j<m;j++) {
            if( str[j] == ‘.‘ )
                num[i][j] = (cnt++);
            else num[i][j] = -1;
        }
    }
    for(int i=0;i<n;i++)
        for(int j=0;j<m;j++) {
            if( num[i][j] == -1 ) continue;
            for(int k=0;k<4;k++) {
                int x = i + dx[k], y = j + dy[k];
                if( x < 0 || y < 0 || x >= n || y >= m ) continue;
                if( num[x][y] == -1 ) continue;
                mat[num[i][j]][num[i][j]]++, mat[num[i][j]][num[x][y]]--;
            }
        }
    printf("%d\n", gauss(cnt - 1, cnt - 1));
}

然后是有向图生成树计数模板题。
bzoj - 4894:天赋/bzoj - 5297: [Cqoi2018]社交网络
两者代码细节上有些不同，这里给出 5297 的代码。

#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN = 250;
const int MOD = 10007;
int pow_mod(int b, int p) {
    int ret = 1;
    while( p ) {
        if( p & 1 ) ret = ret*b%MOD;
        b = b*b%MOD;
        p >>= 1;
    }
    return ret;
}
int mat[MAXN][MAXN];
int gauss(int r, int c) {
    int nwr = 1, nwc = 1, f = 1;
    for(int i=1;i<r;i++)
        for(int j=1;j<c;j++)
            mat[i][j] = (mat[i][j] + MOD)%MOD;
    while( nwr < r && nwc < c ) {
        int mxr = nwr;
        for(int i=nwr+1;i<r;i++)
            if( mat[i][nwc] > mat[mxr][nwc] )
                mxr = i;
        if( mat[mxr][nwc] ) {
            if( nwr != mxr ) {
                for(int i=nwc;i<c;i++)
                    swap(mat[nwr][i], mat[mxr][i]);
                f *= -1;
            }
            int inv = pow_mod(mat[nwr][nwc], MOD-2);
            for(int i=nwr+1;i<r;i++) {
                int x = 1LL*inv*mat[i][nwc]%MOD;
                for(int j=nwc;j<c;j++)
                    mat[i][j] = (mat[i][j] + MOD - 1LL*x*mat[nwr][j]%MOD)%MOD;
            }
            nwr++;
        }
        nwc++;
    }
    int ans = 1;
    for(int i=1;i<r;i++)
        ans = 1LL*ans*mat[i][i]%MOD;
    return (ans*f + MOD)%MOD;
}
int main() {
    int n, m; scanf("%d%d", &n, &m);
    for(int i=1;i<=m;i++) {
        int a, b; scanf("%d%d", &a, &b); a--, b--;
        mat[a][a]++, mat[a][b]--;
    }
    printf("%d\n", gauss(n, n));
}

然后是一道不是特别模板的题。
bzoj - 3534: [Sdoi2014]重建
给出每条边的出现概率，最后出现的图是树的概率。

可以发现这个概率其实等于所有树边出现的概率*所有非树边不出现的概率。
但是因为我们求解生成树只能统计树边。所以我们稍微转换一下，变成：所有树边出现的概率*所有边不出现的概率/所有树边不出现的概率。
然后就可以跑带权的矩阵树定理了。

注意涉及到除法，除数可能为 0（某条边一定为 0）。我们给它微小扰动一下（给一定出现的边的概率减去 eps）。

最后……如果你输出 5 位小数都是过不了的，反正我直接输出 10 位小数。

#include<cmath>
#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN = 50;
const double EPS = 1E-7;
double mat[MAXN][MAXN];
double gauss(int r, int c) {
    int nr = 0, nc = 0; double f = 1;
    while( nr < r && nc < c ) {
        int mxr = nr;
        for(int i=nr+1;i<r;i++)
            if( fabs(mat[i][nc]) > fabs(mat[mxr][nc]) ) mxr = i;
        if( fabs(mat[mxr][nc]) > EPS ) {
            if( nr != mxr ) {
                for(int j=nc;j<c;j++)
                    swap(mat[nr][j], mat[mxr][j]);
                f *= -1;
            }
            for(int i=nr+1;i<r;i++) {
                double x = mat[i][nc]/mat[nr][nc];
                for(int j=nc;j<c;j++)
                    mat[i][j] = (mat[i][j] - x*mat[nr][j]);
            }
            nr++;
        }
        nc++;
    }
    double ans = 1;
    for(int i=0;i<r;i++)
        ans = ans*mat[i][i];
    return ans*f;
}
int main() {
    double f = 1;
    int n; scanf("%d", &n);
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++) {
            double x; scanf("%lf", &x);
            if( x == 1 ) x = 1 - EPS;
            mat[i][i] += x/(1-x), mat[i][j] -= x/(1-x);
            if( i < j ) f *= (1-x);
        }
    printf("%.10lf\n", gauss(n-1, n-1)*f);
}

然后可能不是模板题的一道题。
bzoj - 4596: [Shoi2016]黑暗前的幻想乡
给出 n-1 个公司以及每个公司能够修出来的公路，现在让你统计有多少种修公路的方案，使得每个公司都能负责其中的恰好一条公路且 n 个点连通。
注意如果两种方案即使长得一样，修建某一条公路的公司如果不一样我们也认为是不同的方案。n<= 17。

初看毫无思路。但是如果我们消去每个公司恰好负责一条路这个限制就是一个很简单的生成树计数了。
因此，也许我们可以容斥？
一个方案不合法，则至少有一个公司没有参与到修建公路的工程中。
我们用总方案 - 至少有一个公司没参与 + 至少有两个公司没参与 - ...。
二进制枚举，暴力重构矩阵跑高斯消元。

#include<vector>
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
typedef pair<int, int> pii;
const int MOD = int(1E9) + 7;
const int MAXN = 17;
int pow_mod(int b, int p) {
    int ret = 1;
    while( p ) {
        if( p & 1 ) ret = 1LL*ret*b%MOD;
        b = 1LL*b*b%MOD;
        p >>= 1;
    }
    return ret;
}
int A[MAXN][MAXN];
int gauss(int r, int c) {
    int nwr = 0, nwc = 0, f = 1;
    for(int i=0;i<r;i++)
        for(int j=0;j<c;j++)
            A[i][j] = (A[i][j] + MOD)%MOD;
    while( nwr < r && nwc < c ) {
        int mxr = nwr;
        for(int i=nwr+1;i<r;i++)
            if( A[i][nwc] > A[mxr][nwc] ) mxr = i;
        if( A[mxr][nwc] ) {
            if( mxr != nwr )
                swap(A[mxr], A[nwr]), f *= -1;
            for(int i=nwr+1;i<r;i++)
                while( A[i][nwc] ) {
                    if( A[nwr][nwc] ) {
                        int x = A[i][nwc]/A[nwr][nwc];
                        for(int j=nwc;j<c;j++)
                            A[i][j] = (A[i][j] + MOD - 1LL*x*A[nwr][j]%MOD)%MOD;
                    }
                    swap(A[i], A[nwr]), f *= -1;
                } 
            nwr++;
        }
        nwc++;
    }
    int ans = 1;
    for(int i=0;i<r;i++)
        ans = 1LL*ans*A[i][i]%MOD;
    return ans*f;
}
int bits[1<<MAXN];
vector<pii>edge[MAXN];
int main() {
    int N; scanf("%d", &N);
    for(int i=0;i<N-1;i++) {
        int m, x, y; scanf("%d", &m);
        for(int j=1;j<=m;j++) {
            scanf("%d%d", &x, &y);
            edge[i].push_back(make_pair(x-1, y-1));
        }
    }
    int ans = 0, t = (1<<(N-1));
    for(int s=0;s<t;s++) {
        for(int i=0;i<N;i++)
            for(int j=0;j<N;j++)
                A[i][j] = 0;
        for(int i=0;i<N-1;i++)
            if( (1<<i) & s ) {
                for(int j=0;j<edge[i].size();j++) {
                    int u = edge[i][j].first, v = edge[i][j].second;
                    A[u][u]++, A[v][v]++, A[u][v]--, A[v][u]--;
                }
            }
        bits[s] = bits[s>>1] + (s & 1);
        int f = ( (N - bits[s] - 1) & 1 ) ? -1 : 1;
        ans = (ans + f*gauss(N-1, N-1))%MOD;
    }
    printf("%d\n", (ans + MOD)%MOD);
}

@算法 - 7@ matrix - tree 定理（矩阵树定理）

标签：算法   fabs   put   end   push   有一个   了解   %s   答案   

原文地址：https://www.cnblogs.com/Tiw-Air-OAO/p/10280720.html