P4111 [HEOI2015]小 Z 的房间(矩阵-树定理)
矩阵-树定理
(Kirchhoff's matrix tree theorem)
前言:
首先这个玩意集训的时候就已经讲过了。
但是当时并没有认真听。
现在回头看这玩意……
前置知识:
\(2\):行列式 极其重要
\(3\):高斯消元 极其重要
\(4\):矩阵树定理 这个了解一下这个名字就行了,要不然你直接看会了还看我的博客干屁 不重要
\(5\):柯西–比内公式 这个其实不用掌握
\(6\):一些奇怪的人名:
矩阵-树定理内容:
官方定义:矩阵-树定理是一个计数定理.若连通图 \(G\) 的邻接矩阵为 \(A\) ,将 \(A\) 的对角线 \((i,i)\) 元素依次换为节点 \(i\) 的度 \(d(i)\),其余元素 \((i,j) (j!=i)\) 取 \(A(i,j)\) 的相反数,所得矩阵记为 \(M\) ,则 \(M\) 的每个代数余子式相等,且等于 \(G\) 的生成树的数目.这就是矩阵一树定理.我们常常称矩阵 \(M\) 为基尔霍夫矩阵。
其实这玩意还有一些其他的说法:
比如比较阳间的说法一:将拉普拉斯矩阵去掉任意的一行和一列,得到的矩阵求行列式,即是原图的生成树数量。(拉普拉斯矩阵下面会说到)
还有更通俗易懂的说法二:(度数矩阵(对角线矩阵)-邻接矩阵)的行列式等于该图生成树个数(行列式下面也会说到)。
还有一些奇奇怪怪的说法三:对于一个无向图 \(G\) ,它的生成树个数等于其基尔霍夫 Kirchhoff 矩阵任何一个 \(N−1\) 阶主子式的行列式的绝对值。
其实他们说的是一回事,完全一样。
随便记一个就行了。
拉普拉斯矩阵 (Laplacian Matrix)
其实,它还有一个别的名字:基尔霍夫矩阵(基尔霍夫就是矩阵-树定理的提出者)。
设无向图 \(G\)(允许重边,可以扩展到加权,但是不能有自环) 有 \(n\) 个节点。拉普拉斯矩阵 \(L:=(l_{i,j})_{n\times n}\) 定义如下:
-
若 \(i==j\) 则 \(l_{i,j}=\deg(v_i)\) , \(\deg(v_i)\) 为定点 \(v_i\) 的度。(如果是有向图,那么我们只考虑出度或者入度中的一个).
-
若 $i\neq j $ ,但是顶点 \(v_i\) 和 \(v_j\) 之间有边相连,那么 \(l_{i,j}=-1\).
-
其他情况 \(l_{i,j}=0\).
根据说法一,我们去掉拉普拉斯矩阵的最后一列和最后一行,然后求行列式就得到了问题的解。
说白了,我们首先求出度数矩阵 \(D\)(也就是对角线矩阵,具体意义是这个点有多少个另外的点与他相连,那么对应位置就是几)。
然后求出邻接矩阵 \(A\)(就是表明点与点的连接关系的那个,具体可以参考一些矩阵快速幂的题或者联想一下 floyd)。
之后求出拉普拉斯矩阵(基尔霍夫矩阵) \(W=D-A\).
然后将 \(W\) 去掉任意一行和一列得到 \(K'\) (一般为了方便我们直接去除最后一行和最后一列)。
那么生成树的个数就等于 \(\det(K')\) (就是 \(K'\) 的行列式)。
下面还会提到有向图的拉普拉斯矩阵以及如何推广到到带权。
行列式
一个矩阵对应一个行列式(就是一个数字)。
计算行列式就是需要矩阵化成上三角或下三角,然后观察对角线的元素,如果其中有个一元素为 \(0\) 则整体为 \(0\),否则行列式的值就是对角线上各个元素的乘积。
首先需要知道一些性质:
-
一个行列式和它的转置行列式相等.
-
互换行列式的两行或者两列,行列式的值变号.
-
如果行列式中有两行或者两列完全相同,那么行列式为零.
-
行列式中某一行或者一列中所有元素同乘同一个数 \(k\) ,等于用 \(k\) 乘以行列式.
-
行列式中如果有两行或者或者两列元素对应成比例,那么行列式为零.
-
把行列式的某一行或者某一列的各个元素乘以同一个数以后加到另外一行或者一列的对应元素上去,行列式的值不变.
以上加黑部分提示我们应该使用高斯消元来转化成三角,然后求得行列式。
需要注意的是一定是三角而不是对角线矩阵。
于是特别喜欢 高斯-约旦 消元的同学还是要学一学经典高斯的。
在把矩阵化成三角时需要注意,在题目中模数可能为任意,这样在剩余系就不一定存在逆元,所以需要使用辗转相除法。
或者如果这道题的模数是质数的话也可以使用乘法逆元。
消元结束以后直接算出行列式的值就好了。
int ans=1;
for(re int i=1;i<=tmp;i++){
for(re int j=i+1;j<=tmp;j++){
while(a[j][i]){
int t=a[i][i]/a[j][i];
for(re int k=i;k<=tmp;k++)a[i][k]=(a[i][k]-a[j][k]*t%mod+mod)%mod;
swap(a[i],a[j]);
ans*=-1;
}
}
ans=(ans*a[i][i]%mod+mod)%mod;
}
return (ans%mod+mod)%mod;
这里同时给出对于模数是质数时,求解行列式的方法:
inline int det(){
int ans=1;
for(re int i=1;i<=tmp;i++){
for(re int j=i;j<=tmp;j++){
if(a[j][i]){
swap(a[i],a[j]);
if(i!=j)ans=mod-ans;
break;
}
}
int t=qpow(a[i][i],mod-2);
for(re int j=i+1;j<=tmp;j++){
if(a[j][i]){
for(re int k=tmp;k>=i;k--){
a[j][k]=(a[j][k]-a[i][k]*t%mod*a[j][i]%mod+mod)%mod;
}
}
}
}
for(re int i=1;i<=tmp;i++)ans=ans*a[i][i]%mod;
return ans;
}
需要处理一下负数情况,否则会出问题:
inline void check(){
for(re int i=1;i<=n;i++){
for(re int j=1;j<=n;j++){
if(a[i][j]<0) a[i][j]+=mod;
}
}
}
对于本题,并不是所有的点都可以用来求行列式,只有不是墙的点才能放进矩阵。
否则根据上面的性质:如果其中有个一元素为 \(0\) 则整体为 \(0\) 。
就会直接 WA 了。
对于加边,只加向下或者向右的边,可以保证不重复。
for(re int i=1;i<=n;i++){
for(re int j=1;j<=m;j++){
cin>>a[i][j];
if(a[i][j]=='.')id[i][j]=++cnt;//用cnt来重新编号
}
}
本题相关:
\(1\):这个题需要注意,一定要开 long long.
我们可以观察一下这个题的样例一:
构建出来拉普拉斯矩阵之后并且去掉最后一行和最后一列,矩阵是这样的:
\(\begin{bmatrix} -1&0&2\\0&1&1\\0&0&999999996\end{bmatrix}\)
然后 ,不开 long long 就没有然后了。
\(2\):构建完拉普拉斯矩阵以后千万不要忘记减掉最后一行和最后一列(不要上头。
\(3\):上三角,上三角,上三角!(针对某些高斯-约旦党。
\(4\):取 \(mod\) 的时候注意处理负数。
复杂度分析:
由于我们使用了辗转相除法来进行高斯消元,所以原来 \(O(n^3)\) 的高斯消元被加了一个 \(log\) ,时间复杂度 \(O(n^3m^3log(nm))\)
CODE:
//#define LawrenceSivan
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define re register
const int N=150;
const int mod=1e9;
#define INF 0x3f3f3f3f
#define int long long//记得long long
struct matrix {//这里是矩阵的板子
int n,m;
int tmp;
int a[N][N];
matrix(){}
matrix(const int &_n,const int &_m):n(_n+1),m(_m+1){
memset(a,0,sizeof(a));
}
inline void clear(){
memset(a,0,sizeof(a));
for(int i=0;i<=n;i++){
a[i][i]=1;
}
}
inline void id(){
for(int i=0;i<=n;i++){
a[i][i]=1;
}
}
inline matrix operator + (const matrix &b)const{
matrix res(n,b.m);
for(int i=0;i<=n;i++){
for(int j=0;j<=m;j++){
res.a[i][j]=(a[i][j]+b.a[i][j])%mod;
}
}
return res;
}
inline matrix operator - (const matrix &b)const{
matrix res(n,b.m);
for(int i=0;i<=n;i++){
for(int j=0;j<=m;j++){
res.a[i][j]=(a[i][j]-b.a[i][j])%mod;
}
}
return res;
}
inline matrix operator * (int c)const{
matrix res(n,m);
for(int i=0;i<=n;i++){
for(int j=0;j<=m;j++){
res.a[i][j]=(c*a[i][j])%mod;
}
}
return res;
}
inline matrix operator * (const matrix &b)const{
matrix res(n,b.m);
for(int i=0;i<=n;i++){
for(int j=0;j<=b.m;j++){
for(int k=0;k<=m;k++){
res.a[i][j]=(res.a[i][j]+a[i][k]*b.a[k][j])%mod;
}
}
}
return res;
}
inline matrix operator ^ (int p)const{
matrix res(n,m),x = *this;
res.clear();
for(;p;p>>=1,x=x*x){
if(p&1)res=res*x;
}
return res;
}
inline matrix trans(){
matrix res(m,n);
for(int i=0;i<=n;i++){
for(int j=0;j<=m;j++){
res.a[i][j]=a[j][i];
}
}
return res;
}
inline int det(){//高斯消元与求行列式
int ans=1;
for(re int i=1;i<=tmp;i++){
for(re int j=i+1;j<=tmp;j++){
while(a[j][i]){//辗转相除
int t=a[i][i]/a[j][i];
for(re int k=i;k<=tmp;k++){
a[i][k]=(a[i][k]-a[j][k]*t%mod+mod)%mod;
}
swap(a[i],a[j]);
ans*=-1;//交换以后需要变号
}
}
ans=(ans*a[i][i]%mod+mod)%mod;//勤处理负数
}
return (ans%mod+mod)%mod;
}
}W,D,A;
//W 基尔霍夫矩阵(拉普拉斯矩阵)
//D 度数矩阵(对角线矩阵)
//A 邻接矩阵
int n,m;
char a[15][15];
int id[15][15],cnt;//id 用于重新编号
inline void add(int u,int v){
A.a[u][v]=1;
A.a[v][u]=1;
}
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){x=x*10+(ch^48);ch=getchar();}
return x*f;
}
signed main(){
#ifdef LawrenceSivan
freopen("aa.in","r",stdin);
freopen("aa.out","w",stdout);
#endif
n=read();m=read();
for(re int i=1;i<=n;i++){
for(re int j=1;j<=m;j++){
cin>>a[i][j];
if(a[i][j]=='.')id[i][j]=++cnt;
}
}
W=matrix(cnt,cnt);//初始化矩阵
D=matrix(cnt,cnt);
A=matrix(cnt,cnt);
for(re int i=1;i<=n;i++){//构造邻接矩阵
for(re int j=1;j<=m;j++){
if(id[i][j]&&id[i][j-1])add(id[i][j],id[i][j-1]);
if(id[i][j]&&id[i-1][j])add(id[i][j],id[i-1][j]);
}
}
for(re int i=1;i<=cnt;i++){//构造度数矩阵
for(re int j=1;j<=cnt;j++){
if(A.a[i][j]) D.a[i][i]++;
}
}
W=D-A;//两者作差得到基尔霍夫矩阵(拉普拉斯矩阵)
W.tmp=--cnt;//除掉最后一行与最后一列
/*for(re int i=1;i<=cnt;i++){
for(re int j=1;j<=cnt;j++){
cout<<W.a[i][j]<<" ";
}
cout<<endl;
}*/
printf("%lld\n",W.det());//求行列式
return 0;
}
关于拉普拉斯矩阵的有向图情况:
像我们上面所说的,只考虑入度或者出度就可以了。
但是考虑这两种是不一样的。
设有向图 \(G\)(允许重边,可以扩展到加权,但是不能有自环) 有 \(n\) 个节点。出度矩阵 \(D=(l_{i,j})_{n\times n}\) 定义如下:
-
若 \(i==j\) 则 \(l_{i,j}=\deg^{out}(v_i)\) , \(\deg^{out}(v_i)\) 为定点 \(v_i\) 的出度。
-
若 $i\neq j $ ,则 \(l_{i,j}=0\)
类似地可以定义入度矩阵。
定义 \(edge(i,j)\) 为点 \(i\) 指向点 \(j\) 的有向边数。
定义邻接矩阵 \(A=(l_{i,j})_{n\times n}\)
- 若 \(i\neq j\) , \(l_{i,j}=edge(i,j)\)
于是 就可以定义 \(Laplace\) 矩阵 \(L^{out}\) 和 \(L^{in}\) 。
-
当为 \(L^{in}\) 时,统计的是外向树。(所有边指向叶节点)
-
当为 \(L^{out}\) 时,统计的是内向树。(所有边指向根节点)
内向是出度,外向是入度(一定要记牢,不要混淆)
如果指定了根节点为\(root\) 那么要在计算行列式的时候要去掉的就不是任意一行和一列了,而是第 \(root\) 行和第 \(root\) 列。
inline void add(int t,int u,int v,int w){
if(t==0){//内向树,出度矩阵
(D.a[u][u]+=w)%=mod;
(A.a[u][v]+=w)%=mod;
}
else{//外向树,入度矩阵
(D.a[v][v]+=w)%=mod;
(A.a[u][v]+=w)%=mod;
}
}
因此如果要统计一张图所有的内向树形图,只要枚举所有的根 \(root\) 并对所有的行列式求和即可。
如果要统计一张图所有的外向树形图,只要枚举所有的根 \(root\) 并对所有的行列式求和即可。
关于拉普拉斯矩阵扩展到带权
- 求生成树边权积之和
将上述中的 \(D\) , \(A\) 改为记录相应边边权和,则行列式从生成树个数变为对应的生成树所有边权之积的和。
于是可以发现,只要把边权全都看成 \(1\) ,就是上面这道题我们解决的情况:求生成树个数。