矩阵学习笔记
矩阵是一种数学概念,在 \(OI\) 中有着重要应用。
一个矩阵有行,列,以及里面的数字。如图便是一个 \(2\) 行 \(3\) 列的矩阵:
矩阵数乘
\(\lambda A\) 就是将 \(\lambda\) 依次乘进每个矩阵元素。
矩阵乘法
\(A\times B=C\),那么 \(C_{i,j}=\sum A_{i,k}\times B_{k,j}\)。也就是说,\(n_1,m_1\) 规格的矩阵 \(A\),和 \(n_2,m_2\) 规格的矩阵 \(B\),相乘得到 \(n_2,m_2\) 规格的矩阵 \(C\),当且仅当在 \(m_1=n_2\) 的情况下,两个矩阵可以相乘。
- 单位矩阵:对角线全是 \(1\) 的矩阵 \(E\),称为单位矩阵,满足 \(E\times A=A\)。如图:
- 矩阵乘法满足结合律,但是不满足交换律。
广义矩阵乘法
其中 \(\otimes,\oplus\) 分别表示两种运算符。
若 \(\otimes\) 对 \(\oplus\) 有分配律,那么这种矩阵运算就具有结合律。
比如,最基本的矩阵乘法,\(\otimes=\times,\oplus=+\),有 \(a\times(b+c)=a\times b+a\times c\),有结合律。
图论常用的变形,\(\otimes=+,\oplus=\min\),有 \(a+\min(b,c)=\min(a+b,a+c)\),所以也有结合律。类似的还有很多。
矩阵题目类型
优化DP(通常为线性递推式)
矩阵加速(数列)
显然有如下关系:
由于矩阵具有结合率,所以我们将 \(0,1\) 矩阵做快速幂即可。
[HNOI2011] 数学作业
技巧:分段处理。
首先思考朴素怎么做,有如下式子:
其中 \(k\) 表示数字 \(i\) 的位数。最终的答案就是:\(f_n\)。
可以式子用矩阵表示为:
关键在于,\(k\) 是会变的,但是 \(n\le10^{18}\),于是 \(k\in [1,18]\),那么我们按照 \(k\) 这 \(18\) 个值分别处理即可。
Addition Robot
首先观察到 \((A,B)\rightarrow (A+B,A)\) 这样的式子实际上属于线性递推式,考虑用矩阵乘法的方式进行表达。
所以我们只需要维护一段区间的矩阵乘积即可,由于此处的矩阵并不相同,所以实际上区间\((L,R)\)乘到 \((A,B)\) 的矩阵,是从 \((R,L)\) 乘在一起的矩阵,也就是反过来。
用线段树维护区间信息,我们设 \(tr_{x,0}\) 表示节点 \(x\) 当前的情况,\(tr_{x,1}\) 表示节点 \(x\) 与当前相反的情况。
对于一次区间取翻操作,我们只需将两者交换,打懒标记,并向上更新即可。由于需要维护从右到左的乘积,所以用右儿子乘左儿子即可。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
LL read() {
LL sum=0,flag=1; char c=getchar();
while(c<'0'||c>'9') {if(c=='-') flag=-1; c=getchar();}
while(c>='0'&&c<='9') {sum=sum*10+c-'0'; c=getchar();}
return sum*flag;
}
const LL MOD=1e9+7;
const int N=1e5+10;
int n,q;
string s;
struct Matrix {
int n,m;
LL mx[3][3];
Matrix() {n=m=0; memset(mx,0,sizeof(mx));}
};
Matrix mul(Matrix a,Matrix b) {
Matrix c; c.n=a.n; c.m=b.m;
for(int i=1;i<=a.n;i++) {
for(int j=1;j<=b.m;j++) {
for(int k=1;k<=a.m;k++) {
c.mx[i][j]=(c.mx[i][j]+a.mx[i][k]*b.mx[k][j]%MOD)%MOD;
}
}
}
return c;
}
Matrix init(int n) {
Matrix c; c.n=c.m=n;
for(int i=1;i<=n;i++) c.mx[i][i]=1;
return c;
}
Matrix ksm(Matrix a,LL b) {
Matrix c=init(a.n);
while(b) {
if(b&1) c=mul(c,a);
a=mul(a,a);
b>>=1;
}
return c;
}
Matrix tr[N<<2][2];
int st[N<<2],tag[N<<2];
void build(int nd,int l,int r) {
if(l==r) {
Matrix a; a.n=a.m=2;
a.mx[1][1]=a.mx[1][2]=a.mx[2][2]=1;
Matrix b; b.n=b.m=2;
b.mx[1][1]=b.mx[2][1]=b.mx[2][2]=1;
if(s[l]=='A') {
tr[nd][0]=a;
tr[nd][1]=b;
}
else {
tr[nd][0]=b;
tr[nd][1]=a;
}
return ;
}
int mid=l+r>>1;
build(nd<<1,l,mid); build(nd<<1|1,mid+1,r);
tr[nd][0]=mul(tr[nd<<1|1][0],tr[nd<<1][0]);
tr[nd][1]=mul(tr[nd<<1|1][1],tr[nd<<1][1]);
}
void add(int nd) {
swap(tr[nd][0],tr[nd][1]);
tag[nd]^=1;
}
void pushdown(int nd) {
if(!tag[nd]) return ;
add(nd<<1); add(nd<<1|1);
tag[nd]=0;
}
void change(int nd,int l,int r,int x,int y) {
if(l>y||r<x) return ;
if(l>=x&&r<=y) return add(nd);
pushdown(nd);
int mid=l+r>>1;
change(nd<<1,l,mid,x,y);
change(nd<<1|1,mid+1,r,x,y);
tr[nd][0]=mul(tr[nd<<1|1][0],tr[nd<<1][0]);
tr[nd][1]=mul(tr[nd<<1|1][1],tr[nd<<1][1]);
}
Matrix query(int nd,int l,int r,int x,int y) {
if(l>y||r<x) return init(2);
if(l>=x&&r<=y) return tr[nd][0];
pushdown(nd);
int mid=l+r>>1;
return mul(query(nd<<1|1,mid+1,r,x,y),query(nd<<1,l,mid,x,y));
}
int main() {
cin>>n>>q>>s;
s=" "+s;
build(1,1,n);
while(q--) {
int opt,l,r,x,y;
cin>>opt>>l>>r;
if(opt==1) {
change(1,1,n,l,r);
}
else {
cin>>x>>y;
Matrix a,b=query(1,1,n,l,r);
a.n=2; a.m=1; a.mx[1][1]=x; a.mx[2][1]=y;
b=mul(b,a);
cout<<b.mx[1][1]<<" "<<b.mx[2][1]<<'\n';
}
}
return 0;
}
矩阵本身的计算
[Cnoi2021] 矩阵
由于 \(n\le 10^5\),所以甚至不能将整个矩阵存下来,更不能进行计算。
由于矩阵 \(A_{i,j}=a_i\times b_j\),所以即可转化成如下矩阵形式:
这样乘出来的矩阵大小是 \(n\times n\) 的,但是若用 \(b\times a\),则矩阵大小便是 \(1\)。所以我们可以利用这个特点转化:
对中间的做普通快速幂即可。
这种计算方法也常常用于简化,常见方式:用 \(n\times n\) 的矩阵与 \(n\times 1\) 的矩阵相乘,得到的仍是 \(n\times 1\) 的矩阵,并且单次时间复杂度为 \(\rm O(n^2)\)。
与图论的结合
这类问题通常带有恰好操作 \(k\) 次的字眼。
[USACO07NOV] Cow Relays G
如果我们求出来一个恰好经过 \(x\) 条边的矩阵 \(f1\),其中 \(f1_{i,j}\) 表示 \(i\sim j\) 的的最短路,以及恰好经过 \(y\) 条边的矩阵 \(f2\),按照以下式子计算出 \(f3\)。
for(int i=1;i<=a.n;i++) {
for(int j=1;j<=b.m;j++) {
for(int k=1;k<=a.m;k++) {
c.mx[i][j]=min(c.mx[i][j],a.mx[i][k]+b.mx[k][j]);
}
}
}
那么 \(f3\) 计算结果便是恰好经过 \(x+y\) 条边的矩阵。
同时也论证过 \(\min\) 和 \(+\) 之间的广义矩阵乘法是满足结合率的,所以我们设置恰好经过一条边的矩阵,然后按照上面的代码做一遍快速幂即可。
类似问题:[SCOI2009] 迷路
[SDOI2009] HH去散步
本题需要满足边与边之间的关系,并不是点与点之间的关系,所以我们为了方便处理,将边转化为点,进行考察。
对于这张图,我们将边标号之后,其中 \(1,2\) 实际为一条双向边,那么即可在 \(1,3\) 以及 \(2,4\) 之间连边,为了满足题目要求,\(1,2\) 以及 \(3,4\) 之间不连边。
那么先预处理出矩阵后,做矩阵快速幂即可。
[NOI Online #3 提高组] 魔法值
异或的形式看上去十分丑陋,我们先想办法进行转化。
首先可以想到,为了方便计算,肯定是将所有点都加进 \(f_{x,j}\) 的计算式当中,那么对于与 \(x\) 号城市不相连的点,我们显然可以采用 \(\times 0\) 的方式将其干掉,那么从广义矩乘的角度来看,便是里面 \(\times\),外面 \(\oplus\),但是显然,\(\times\) 对 \(\oplus\) 不符合分配律,也就是说, \(a\times (b\oplus c)=(a\times b)\oplus (a\times c)\)。
但是,乘法的对象只有 \(0\) 和 \(1\),不难证明,在这种情况下,是满足分配律的。
所以我们便可预处理出 \(0,1\) 矩阵,然后快速幂,单词询问时间复杂度 \(\rm O(n^3\log a)\),总时间复杂度即为 \(\rm O(qn^3\log a)\),无法通过本题。
考虑优化。
注意到 \(f\) 矩阵实际上只有 \(n\times 1\),关键在于 \(01\) 矩阵(记做 \(A\))是 \(n\times n\),且转移的过程用 \(A\times f\),所以我们可以采用倍增预处理,先预处理出 \(A^{2^k}\),时间复杂度 \(\rm O(n^3\log a)\),并且每次处理询问,用预处理的矩阵乘 \(f\),乘出来的仍是 \(n\times 1\),单次询问时间复杂度 \(\rm O(n^2\log a)\),总时间复杂度为 \(\rm O(qn^2\log a+n^3\log a))\),可以通过本题。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
LL read() {
LL sum=0,flag=1; char c=getchar();
while(c<'0'||c>'9') {if(c=='-') flag=-1; c=getchar();}
while(c>='0'&&c<='9') {sum=sum*10+c-'0'; c=getchar();}
return sum*flag;
}
const int N=1e5+10;
const LL S=(1ll<<32)-1;
int n,m,q;
struct node {
LL a; int id;
}s[N];
LL ans[N];
struct Matrix {
int n,m;
LL mx[110][110];
Matrix() {n=m=0; memset(mx,0,sizeof(mx));}
}bz[35];
Matrix mul(Matrix a,Matrix b) {
Matrix c; c.n=a.n; c.m=b.m;
for(int i=1;i<=a.n;i++) {
for(int j=1;j<=b.m;j++) {
for(int k=1;k<=a.m;k++) {
c.mx[i][j]^=a.mx[i][k]&b.mx[k][j];
}
}
}
return c;
}
Matrix init(int n) {
Matrix c; c.n=c.m=n;
for(int i=1;i<=n;i++) c.mx[i][i]=1;
return c;
}
Matrix ksm(Matrix a,LL b) {
Matrix c;
while(b) {
if(b&1) {
if(!c.m) c=a;
else c=mul(c,a);
}
a=mul(a,a);
b>>=1;
}
return c;
}
int cmp1(node x,node y) {
return x.a<y.a;
}
int main() {
// freopen("a.in","r",stdin);
// freopen("a.out","w",stdout);
cin>>n>>m>>q;
Matrix f; f.n=n; f.m=1;
for(int i=1;i<=n;i++) cin>>f.mx[i][1];
Matrix zy; zy.n=zy.m=n;
bz[0].m=bz[0].n=n;
for(int i=1;i<=m;i++) {
int u=read(),v=read();
bz[0].mx[u][v]=bz[0].mx[v][u]=S;
}
for(int i=1;i<=32;i++) bz[i]=mul(bz[i-1],bz[i-1]);
while(q--) {
LL x; cin>>x;
Matrix ans=f;
for(int i=32;i>=0;i--) {
if(x>=(1ll<<i)) {
ans=mul(bz[i],ans);
x-=(1ll<<i);
}
}
cout<<ans.mx[1][1]<<endl;
}
return 0;
}