浅谈矩阵在OI中的应用
本文作者太菜了，不会用专业的数学术语表达，在本文中用相对通俗的语言进行讲解。
定义
矩阵（Matrix）是一个按照长方阵列排列的复数或实数集合
可以想象为一个 \(n\times m\) 的二维数组，可以表示为：
如果 \(n=m\)，那么 \(A\) 也被称为 \(n\) 阶矩阵。
单位矩阵 \(I\) 表示 \(I_{i,i}=1\)，其余都为 \(0\) 的矩阵。
基本运算
加减法表示两个大小相同的矩阵每个对于位置相加。
即
一个矩阵乘上一个数，表示每个位置都乘上这个数。
即
乘法
两个矩阵相乘仅当第一个矩阵 \(A\) 的列数和另一个矩阵 \(B\) 的行数相等时才能乘。如 \(A\) 是 \(m\times n\)矩阵和 \(B\) 是 \(n\times p\) 矩阵，它们的乘积 \(C\) 是一个 \(m\times p\) 矩阵。
乘法满足
计该运算为 \(C=AB\)。
也就是 \(C_{i,j}\) 表示 \(A\) 的第 \(i\) 行，\(B\) 的第 \(j\) 列的每一个元素乘起来的和。
Matr operator *(const Matr &x,const Matr &y){
Matr z;
memset(z.cnt,0,sizeof(z.cnt));
for(int k=1;k<=n;++k)
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
z.cnt[i][j]=(z.cnt[i][j]+x.cnt[i][k]*1ll*y.cnt[k][j]%mod)%mod;
return z;
}
需要注意的是，矩阵乘法满足结合律和分配率，但不满足交换率。
因此，矩阵乘法是可以用快速幂的。
需要注意的是，对于任意矩阵 \(A\)，\(A^0=I\)，所以可以将 \(I\) 定义为初值。
板子题，直接套一个矩阵乘法和快速幂就好了。
时间复杂度 \(O(n^3\log k)\)
#include<bits/stdc++.h>
using namespace std;
const long long N=105;
long long n,k,mod=1E9+7;
struct hzy{long long cnt[N][N];}a;
hzy operator *(const hzy &x,const hzy &y){
hzy z;
memset(z.cnt,0,sizeof(z.cnt));
for(long long k=1;k<=n;++k)
for(long long i=1;i<=n;++i)
for(long long j=1;j<=n;++j)
z.cnt[i][j]=(z.cnt[i][j]+x.cnt[i][k]*y.cnt[k][j]%mod)%mod;
return z;
}
hzy kuai(hzy a,long long b){
hzy c;
for(long long i=1;i<=n;i++)c.cnt[i][i]=1;
for(;b;b>>=1){
if(b&1)c=c*a;
a=a*a;
}
return c;
}
int main(){
scanf("%lld%lld",&n,&k);
for(long long i=1;i<=n;i++)
for(long long j=1;j<=n;j++)
scanf("%lld",&a.cnt[i][j]);
a=kuai(a,k);
for(long long i=1;i<=n;i++){
for(long long j=1;j<=n;j++)printf("%lld ",a.cnt[i][j]);
puts("");
}
}//之前写的代码，写的实在太丑了
矩阵快速幂有啥用？
矩阵快速幂的一个重要应用就是优化 dp 转移。
来个例子斐波那契数列。
众所周知，斐波那契数列是可以线性递推的，时间复杂度 \(O(n)\)。
尝试使用矩阵乘法进行优化。
我们发现，\(fib_i\) 只与 \(fib_{i-1}\) 和 \(f_{i-2}\) 有关。
那么我们构造一个矩阵\(\begin{bmatrix} fib_{i-2} & fib_{i-1} \\ \end{bmatrix}\)。
我们通过乘法将它变成\(\begin{bmatrix} fib_{i-1} & fib_{i-2}+fib_{i-1} \\ \end{bmatrix}\)
容易得到它乘上的矩阵\(A=\begin{bmatrix} 0 & 1 \\ 1 & 1 \\ \end{bmatrix}\)
事实上，我们发现 \(A_{i,j}\) 就表示在上个状态中的 \(f_{t-1,i}\) 给 \(f_{t,j}\) 的贡献就是 \(A_{i,j}\times f_{t-1,i}\)，因此很容易推出。
乘一次，就向后转移一次，根据结合律对后面相同的矩阵进行合并。
容易发现我们要运算的实际上是 \(\begin{bmatrix} 1 & 1 \\ \end{bmatrix}\times A^{n-1}\)，直接矩阵快速幂即可。
时间复杂度 \(O(\log n)\)。
mat a,b;
int main(){
scanf("%lld",&n);
a.cnt[1][1]=0;a.cnt[1][2]=1;
a.cnt[2][1]=1;a.cnt[2][2]=1;
a=kuai(a,n-1);
b.cnt[1][1]=1;
b.cnt[1][2]=1;
b=b*a;
printf("%lld",b.cnt[1][1]%mod);
}
关于斐波那契数列的计算，还可以用生成函数优化为 \(O(1)\)。
练习：【模板】矩阵加速。
还是类似，考虑原地不动就是一个自环，自爆就见一个虚点，所有的点都连向它。
设 \(f_{i,j}\) 表示 \(i\) 时刻最后停留在 \(j\) 的方案数。
暴力 dp 过不了，考虑每次转移的本质都相同，直接矩阵快速幂。
转移矩阵实际上就是邻接矩阵（根据上面推出转移矩阵的结论）。
然后直接做就好了。
时间复杂度 \(O(n^3\log t)\)。
int main(){
Matr c;
memset(c.cnt,0,sizeof(c.cnt));
int ans=0;
scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++){
int x,y;
scanf("%d%d",&x,&y);
c.cnt[x][y]=c.cnt[y][x]=1;
}
scanf("%d",&k);
n++;
for(int i=1;i<=n;i++)c.cnt[i][i]=1;
for(int i=1;i<n;i++)c.cnt[i][n]=1;
c=kuai(c,k);
for(int i=1;i<=n;i++)(ans+=c.cnt[1][i])%=mod;
printf("%d",ans);
}
仍然考虑上面的这个问题，如果有 \(Q\) 次询问，每次询问时刻 \(t\) 时的答案怎么做，\(Q\le 100\)。
如果暴力每次做，时间复杂度 \(O(Qn^3\log t)\)，炸了。
容易发现做一次矩阵乘法的时间复杂度是 \(O(n^3)\) 的，但我们发现答案的矩阵是 \(1\times n\) 的。
答案的这个矩阵乘上转移矩阵是 \(O(n^2)\) 的。
设初始矩阵为 \(P\)。
答案就是 \(P\times A^t\)。
快速幂拆成 \(P\times A^{2^{t_0}}\times A^{2^{t_1}}\times A^{2^{t_ 2}} ...\times A^{2^{t_k}}\)
根据结合律，变成\((...(((P\times A^{2^{t_0}})\times A^{2^{t_1}})\times A^{2^{t_ 2}}) ...\times A^{2^{t_k}})\)
每次运算都是 \(O(n^2)\) 的，开始时先预处理后面要乘的 \(A^{2^t}\)。
时间复杂度 \(O(Qn^2\log t+n^3\log t)\)。
练习：
都是推出暴力式子后直接矩阵快速幂的。
定长最短路
考虑矩阵乘法：
Matr operator *(const Matr &x,const Matr &y){
Matr z;
memset(z.cnt,0,sizeof(z.cnt));
for(int k=1;k<=n;++k)
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
z.cnt[i][j]=(z.cnt[i][j]+x.cnt[i][k]*1ll*y.cnt[k][j]%mod)%mod;
return z;
}
你会发现这个东西和 Floyd 不仅形似而且神似。
改亿改，变成 Floyd 的形式：
Matr operator *(const Matr &x,const Matr &y){
Matr z;
memset(z.cnt,0,sizeof(z.cnt));
for(int k=1;k<=n;++k)
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
z.cnt[i][j]=min(z.cnt[i][j],max(x.cnt[i][k]+y.cnt[k][j]));
return z;
}
现在变成了我们熟悉的形式，但有什么用呢？
考虑设邻接矩阵为 \(A\)，那么根据前面 dp 的思路， \(A^L\) 表示走 \(L\) 步后得到的最短路的邻接矩阵。
和上面说的一模一样，需要注意的是邻接矩阵中 \(A_{i,i}=Inf\)，不然就可以一直待在原地。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1005;
struct hzy{int cnt[N][N];}a;
int k,m,st,en,b[N],x[N],y[N],z[N],id[N],tot,dis[N][N];
hzy operator*(const hzy x,const hzy y){
hzy z;
memset(z.cnt,63,sizeof(z.cnt));
for(int k=1;k<=tot;k++)
for(int i=1;i<=tot;i++)
for(int j=1;j<=tot;j++)z.cnt[i][j]=min(z.cnt[i][j],x.cnt[i][k]+y.cnt[k][j]);
return z;
}
hzy kuai(hzy a,int b){
hzy ans;
memset(ans.cnt,63,sizeof(ans.cnt));
for(int i=1;i<=tot;i++)ans.cnt[i][i]=0;
for(;b;b>>=1,a=a*a)if(b&1)ans=ans*a;
return ans;
}
int main() {
memset(dis,63,sizeof(dis));
scanf("%d%d%d%d",&k,&m,&st,&en);
for(int i=1;i<=m;i++){
scanf("%d%d%d",&z[i],&x[i],&y[i]);
b[x[i]]=1,b[y[i]]=1;
}
for(int i=1;i<=1000;i++)if(b[i])id[i]=++tot;
for(int i=1;i<=m;i++)x[i]=id[x[i]],y[i]=id[y[i]];
for(int i=1;i<=m;i++)
dis[y[i]][x[i]]=dis[x[i]][y[i]]=min(dis[x[i]][y[i]],z[i]);
for(int i=1;i<=tot;i++)for(int j=1;j<=tot;j++)a.cnt[i][j]=dis[i][j];
printf("%d",kuai(a,k).cnt[id[en]][id[st]]);
return 0;
}
如果修改题目条件，变成走至多 \(L\) 步的最短路。
直接每个点连一个边权为 \(0\) 的自环，那么就相当于可以一直留在原地，就做完了。
再来个例题Good Vertices。
和前面几乎一样，还是跑 Flody，只不过修改一下就好了,\(A_{i,j}\) 表示 \(i\) 走 \(L\) 步到 \(j\) 的最小时刻。
mat operator*(mat a,mat b){
mat ans;
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)ans.cnt[i][j]=1e9;
for(int k=1;k<=n;k++)
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)ans.cnt[i][j]=min(ans.cnt[i][j],max(a.cnt[i][k],b.cnt[k][j]));
return ans;
}
\(n\) 很小，依然是 Flody，然后考虑一个转移矩阵表示用一次魔法后的邻接矩阵。
考虑使用一次魔法的改变即可。
转移矩阵：
for(int k=1;k<=n;k++)
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
dis[i][j]=min(dis[i][j],dis[i][k]+dis[k][j]);
for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)p.cnt[i][j]=dis[i][j];
for(int i=1;i<=m;i++)//考虑每一条边，把他变成负数
for(int l=1;l<=n;l++)
for(int r=1;r<=n;r++)
p.cnt[l][r]=min(p.cnt[l][r],dis[l][x[i]]+dis[y[i]][r]-z[i]);
练习：[NOI2020] 美食家。
维护数列
矩阵乘法还可以和数据结构一起食用。
来个题Fibonacci-ish II。
考虑没有修改操作，那么直接用莫队，这样就只需要考虑加一个或者删掉一个数就好了。
考虑新加入一个数会让后面的数都向后平移，即 \(a_x\times fib_y\) 变成了 \(a_x\times fib_{y+1}\)。
那么这个变化就相当于乘上了一个 \(fib\) 的矩阵。
那么这个东西就可以直接用线段树维护了。
code。
没想到吧，人傻常数大，这个代码 Time limit exceeded on test 26 了。
考虑优化。
其实不需要延迟标记，直接考虑单点修改，只要可以维护两个区间直接的合并就好了。
原来两个区间的和分别是
合并成
实际上就是右边的乘上 \(Fib^k\)，因为 \(n\) 很小，直接预处理就好了。
code。
动态dp
动态动态规划。
还在学，快了。
To be continue。
行列式
定义一个行列式运算 \(|A|\)，也可以表示为 \(\det(A)\)，\(A\) 是一个 \(n\) 阶矩阵。
\(P\) 表示一个排列，\(S\) 表示全排列集合，\(f(P)\) 表示 \(P\) 的逆序对数量。
也就是每行每列选一个数乘起来，再乘一个容斥系数。
暴力求的时间复杂度为 \(O(n!\times n)\)
来一些性质：
-
如果 \(A\) 是三角矩阵，那么 \(\det(A)=\prod_{i=1}^n A_{i,i}\)，这个比较显然。
-
交换 \(A\) 的两行或者两列，那么 \(\det(A)\) 取相反数。
证明：考虑两行，取得位置为 \(P_{i},P_j\)。
交换之后实际上就是逆序对中交换了两个数。
讨论一下就可以发现逆序对个数的变化一定是个奇数。
所以行列式的值取相反数。
交换两列也差不多。
- 如果 \(A\) 有两行一模一样，那么 \(|A|=0\)。
证明：交换这两行，\(|A|=-|A|\)。
- 基本变换，将矩阵一行上所有的数乘 \(t\)，加到另一行上，\(|A|\) 不变。
这个也比较简单，加过去，通过乘法分配律拆成 \(t+1\) 个矩阵，然后发现后面的矩阵都有两行相等，都是 \(0\) ，所以 \(|A|\) 不变。
发现没有，这个变换和高斯消元一模一样。
可以用高斯消元消为一个上三角，把对角线乘起来就好了。
时间复杂度 \(O(n^3)\)。
#define For(x,l,r) for(int x=l;x<=r;x++)
int det(mat a){
int ans=1,ff=1;
For(i,1,n){
int k=i;
while(!a.cnt[k][i]&&k<=n)k++;
if(k!=i){ff*=-1;For(j,i,n)swap(a.cnt[i][j],a.cnt[k][j]);}//交换时记得取相反数
if(!a.cnt[i][i])return 0;
ans=1ll*ans*a.cnt[i][i]%mod;
int Inv=kuai(a.cnt[i][i],mod-2);
For(j,i+1,n){//普通高斯消元
if(!a.cnt[j][i])continue;
int t=1ll*Inv*a.cnt[j][i]%mod;
For(p,i,n)a.cnt[j][p]=(a.cnt[j][p]-1ll*t*a.cnt[i][p]%mod+mod)%mod;
}
}
return (ans*ff+mod)%mod;
}
然后发现这并不能过【模板】行列式求值。
当 \(p\) 非质数时，原来的做法会有问题。
考虑中间那个消元过程，就是一个不断减掉的过程，可以辗转相减。
#define For(i,l,r) for(int i=l;i<=r;i++)
int det(mat a){
int ans=1,ff=1;
For(i,1,n)For(j,i+1,n){
while(a.cnt[i][i]){
int div=a.cnt[j][i]/a.cnt[i][i];
For(k,i,n)a.cnt[j][k]=(a.cnt[j][k]-1ll*div*a.cnt[i][k]%mod+mod)%mod;
ff=-ff,swap(a.cnt[i],a.cnt[j]);
}
ff=-ff,swap(a.cnt[i],a.cnt[j]);
}
for(int i=1;i<=n;i++)ans=1ll*ans*a.cnt[i][i]%mod;
return (ans*ff+mod)%mod;
}
时间复杂度仍然是 \(O(n^3)\)。
- 范德蒙德行列式
考虑证明，使用归纳法，对于 \(n=2,3\)，可以直接暴力展开证明。
如果 \(n-1\) 成立，考虑对于 \(n\) 阶矩阵，把第 \(i-1\) 行每个数乘上 \(x_1\)，然后让第 \(i\) 行去减掉，倒着做一遍就可以得到：
LGV引理
先来个模板：[NOI2021] 路径交点
这个引理主要解决有向无环图中不相交路径求和问题。
已知一张有向无环图。
\(f_{i,j}\) 表示 \(i\) 走到 \(j\) 的方案数。
杨氏矩阵
还在学，快了。
To be continue。
