行列式、矩阵树定理-526互联

推荐阅读：矩阵树定理(+行列式) - command_block 的博客。

行列式

定义

这个东西一般用于求解图的生成树个数（矩阵树定理）。

称一个大小为 \(n\times n\) 的矩阵 \(A\) 的行列式为 \(\det(A)\)（或 \(|A|\)）。

\[\det(A)=\sum_{p\texttt{是一个大小为n排列}}(-1)^{F(p)}\prod_{i=1}^{n}A[i][p_i] \]

其中 \(F(p)\) 表示排列 \(p\) 的逆序对数。\(\prod_{i=1}^{n}A[i][p_i]\) 可以理解为每一行选一个，每次选的列不与其他行重复，将选中的数乘起来。

性质

单位矩阵 \(I\) 的行列式为 \(1\)，上三角、下三角矩阵的行列式为对角线的值的乘积。
交换矩阵的任意两行，行列式变号。（交换一个排列的两个值，逆序对数的奇偶性一定发生变化）。
如果矩阵的某一行均为 \(0\) 则行列式为值为 \(0\)。（每行选且仅选一个）
若矩阵 \(A\) 的第 \(i\) 行和第 \(j\) 行相同，行列式的值为零。（交换第 \(i\)、\(j\) 行后行列式变号，但实际上矩阵 \(A\) 没有发生变化即行列式不变号，所以行列式值为零）。
若矩阵 \(A\) 的第 \(i\) 行乘以 \(k\)，行列式的值也乘以 \(k\)。（每一行都会选一个）
\( \det(\begin{matrix}a+e &b+f\\c &d\\\end{matrix})= \det(\begin{matrix}a &b\\c &d\\\end{matrix})+ \det(\begin{matrix}e &f\\c &d\\\end{matrix}) \)（结合行列式定义自行理解）
\( \det(\begin{matrix}a+kc &b+kd\\c &d\\\end{matrix})=\det(\begin{matrix}a &b\\c &d\\\end{matrix}) \)（\(\det(\begin{matrix}kc &kd\\c &d\\\end{matrix})=\det(\begin{matrix}c &d\\c &d\\\end{matrix})=0\)）

高斯消元求解行列式

根据性质 \(6\)，我们可以用高斯消元求解行列式，将原矩阵转化为上三角矩阵、顺便记录交换行的次数的奇偶性即可。

至于如何用第 \(i\) 行将 \(a_{j,i}\) 变为零（\(n\ge j>i\)）。可以用辗转相除法。

这样既可以求出任意一个矩阵在模 \(p\) 下的行列式的值了（\(p\) 可以不为质数）。

Code:

typedef long long LL;
#define arr2D vector<vector<int>>
#define arr vector<int>
int solve(arr2D &a,int n,int mo){
    bool Swap=false; LL ans=1;
    for(int i=1;i<=n;i++){
        int tmp=i;
        for(int j=i;j<=n;j++)if(a[j][i]){tmp=j;break;}
        if(!a[tmp][i]){return 0;}
        for(int j=tmp+1;j<=n;j++)if(a[j][i]&&a[j][i]<tmp)tmp=j;
        if(i!=tmp){a[i].swap(a[tmp]);Swap^=1;}
        for(int j=i+1;j<=n;j++){
            if(a[j][i]>a[i][i]){a[i].swap(a[j]),Swap^=1;}
            while(a[j][i]){
                int d=a[i][i]/a[j][i];
                for(int k=i;k<=n;k++)a[i][k]=(a[i][k]+(LL)(mo-d)*a[j][k]%mo)%mo;
                a[i].swap(a[j]),Swap^=1;
            }
        }
        ans=(LL)1LL*ans*a[i][i]%mo;
    }
    if(Swap)ans=(mo-ans)%mo;
    return ans;
}

矩阵树定理

定义

一般用于求解图的生成树问题。（若无特殊说明，默认图为无向图，图中一定没有自环）

度数矩阵 \(D\) 满足 \(d_{i,i}=du_i\) 其他的 \(d_{i,j}=0,i\not=j\)（在无向图中 \(du_i\) 表示与 \(i\) 相连的边数），邻接矩阵 \(A\) 表示连边情况。其基尔霍夫矩阵为 \(K=D-A\)。则该图的生成树个数为矩阵 \(K\) 删去第 \(r\) 行、\(r\) 列后的矩阵的行列式的值。

如果矩阵树定理只能求生成树个数就太废了，实际上矩阵树定理求的所有生成树的边权乘积之和（求生成树个数时所有边权均为 \(1\)）。

因此将度数矩阵 \(D\) 中的 \(d_{i,i}\) 改为与 \(i\) 相连的边权之和，邻接矩阵 \(A\) 中的 \(a_{i,j}\) 改为 \(i\) 和 \(j\) 之间的边权之和。其基尔霍夫矩阵为 \(K=D-A\)。则该图的所有生成树的边权乘积之和为矩阵 \(K\) 删去第 \(r\) 行、\(r\) 列后的矩阵的行列式的值。

因此在用矩阵树定理求解生成树个数时可以支持重边（边 \((i,j)\) 的权值为 \((i,j)\) 在原图中出现的次数）。

有向图的矩阵树定理

邻接矩阵 \(A\) 为有向图的连边情况。

度数矩阵要分类讨论：

度数矩阵 \(D\) 的 \(d_{i,i}=\sum{j=1}^{n}A_{i,j}\) 时求的是外向树生成树个数。
度数矩阵 \(D\) 的 \(d_{i,i}=\sum{j=1}^{n}A_{j,i}\) 时求的是内向树生成树个数。

若以 \(k\) 为根则其生成树个数为基尔霍夫矩阵删去第 \(k\) 行、\(k\) 列的矩阵的行列式的值。

例题

P4336 [SHOI2016] 黑暗前的幻想乡

直接容斥计算即可，\(\sum_{i}^{}(-1)^{n-P(i)}F(i)\)，\(P(i)\) 表示集合 \(i\) 的元素个数，\(F(i)\) 表示包含集合 \(i\) 内的所有建筑商可以修建的边的图的生成树个数。

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef vector<int> arr;
typedef vector<arr> arr2D;
constexpr int mo=1e9+7;

int sol(arr2D a,int n){
    LL ans=1;bool Swap=false;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            a[i][j]%=mo;
    for(int i=1;i<=n;i++){
        int tmp=i;
        for(int j=i;j<=n;j++)if(a[j][i]){tmp=j;break;}
        if(!a[tmp][i])return 0;
        for(int j=tmp+1;j<=n;j++)if(a[j][i]&&a[j][i]<a[tmp][i])tmp=j;
        if(i!=tmp)a[i].swap(a[tmp]),Swap^=1;
        for(int j=i+1;j<=n;j++){
            if(a[j][i]<a[i][i])a[i].swap(a[j]),Swap^=1;
            while(a[j][i]){
                int d=a[i][i]/a[j][i];
                for(int k=i;k<=n;k++)
                    a[i][k]=(a[i][k]+(LL)1LL*(mo-d)*a[j][k])%mo;
                a[i].swap(a[j]),Swap^=1;
            }
        }
        ans=(LL)ans*a[i][i]%mo;
    }
    if(Swap)ans=(mo-ans)%mo;
    return ans;
}

arr2D Val;
vector<pair<int,int>>e[55];
int n;

int Sol(int S){
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            Val[i][j]=0;
    int t=0;
    for(int i=0;i<n-1;i++)if((S>>i)&1){
        t++;
        for(auto j:e[i])
            Val[j.first][j.first]++,
            Val[j.second][j.second]++,
            Val[j.first][j.second]--,
            Val[j.second][j.first]--;
    }
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            Val[i][j]=(mo+Val[i][j]%mo)%mo;
    int ans=sol(Val,n-1);
    if((t&1)!=((n-1)&1))ans=(mo-ans)%mo;
    return ans;
}

int main(){
    cin>>n;
    Val.resize(n+1,arr(n+1));
    for(int i=0;i<n-1;i++){
        int m;cin>>m;e[i].resize(m);
        for(auto& j:e[i])
            cin>>j.first>>j.second;
    }
    int ans=0;
    for(int i=1;i<(1<<n-1);i++)
        (ans+=Sol(i))%=mo;
    cout<<ans<<'\n';
    return 0;
}

P3317 [SDOI2014] 重建

边权改为 \(\frac{p_{i,j}}{1-p_{i,j}}\)，注意 \(p_{i,j}=1\) 的情况，手动加个误差即可。

#include <bits/stdc++.h>
using namespace std;
typedef vector<double> arr;
typedef vector<arr> arr2D;
constexpr double eps=1e-10;
double sol(arr2D &a,int n){
    double ans=1.;
    for(int i=1;i<=n;i++){
        for(int j=i+1;j<=n;j++)if(fabs(a[j][i])>eps){
            double d=a[j][i]/a[i][i];
            for(int k=i;k<=n;k++)a[j][k]-=d*a[i][k];
        }
        ans*=a[i][i];
    }
    return ans;
}
arr2D a;
int n;
int main(){
    cin>>n;
    a.resize(n+1,arr(n+1,0.));
    double pi=1.;
    for(int i=1;i<=n;i++){
        double sum=0;
        for(int j=1;j<=n;j++){
            cin>>a[i][j];a[i][j]+=eps;
            if(i<j)pi*=(1.-a[i][j]);
            if(i^j)a[i][j]=a[i][j]/(1.-a[i][j]);
            sum+=a[i][j],a[i][j]=-a[i][j];
        }
        a[i][i]=sum;
    }
    printf("%.6lf\n",sol(a,n-1)*pi);
    return 0;
}