浅谈矩阵快速幂

pidan123

在谈矩阵快速幂前，我们先说说什么是矩阵。

关于矩阵

顾名思义，矩阵就是把一些数填入一个矩形形式如下：

\[\begin{bmatrix}1&4&2&3&5\\9&6&4&1&7\\3&4&5&6&7\end{bmatrix} \]

所以呢？

所以我们再来介绍一下矩阵的乘法。

矩阵乘法

矩阵乘法就是将两个矩阵相乘，可是我们要怎么运算呢？百度搜了搜从书上得知，矩阵乘法有如下定义：

设A、B是两个矩阵，令C=A×B，那么：

  1. A的列数必须和B的行数相等；

  2. 设A是n×r的矩阵，B是r×m的矩阵；那么A和B的乘积C是一个n×m的矩阵；

  3. \(C_{i,j}=A_{i,1}×B_{1,j}+A_{i,2}×B_{2,j}+A_{i,3}×B_{3,j}+…+A_{i,n}×B_{n,j}=\sum^{n}_{k=1}A_{i,k}×B_{k,j}\)

例如：

\[\begin{bmatrix}1&4\\2&5\\3&6\end{bmatrix}\times{\begin{bmatrix}1&2&3\\4&5&6\end{bmatrix}}=\begin{bmatrix}1\times1+4\times4&1\times2+4\times5&1\times3+4\times6\\2\times1+5\times4&2\times2+5\times5&2\times3+5\times6\\3\times1+6\times4&3\times2+6\times5&3\times3+6\times6\end{bmatrix}=\begin{bmatrix}17&22&27\\22&29&36\\27&36&45\end{bmatrix} \]

由此我们可以写出代码：

struct matrix{
	int n,m,g[N][N];//n为行数,m为列数,g是矩阵,N为矩阵大小

	matrix operator*(const matrix&b)const{
		matrix c;
		memset(c.g,0,sizeof(c.g));//初始化
		for(int i=1;i<=n;i++)
			for(int j=1;j<=m;j++)
				for(int k=1;k<=b.m;k++)
					c.g[i][j]=(c.g[i][j]+g[i][k]*b.g[k][j]%mod)%mod;
		return c;
	}
}

matrix multiply(matrix a,matrix b){
	matrix c;
	memset(c.g,0,sizeof(c.g));
	for(int i=1;i<=a.n;i++)
		for(int j=1;j<=a.m;j++)
			for(int k=1;k<=b.m;k++)
				c.g[i][j]=(c.g[i][j]+a.g[i][k]*b.g[k][j]%mod)%mod;
	return c;
}

那么矩阵乘法满不满足运算律呢？

有的，矩阵乘法满足结合律，证明如下：

设n阶矩阵为：

\(A=(a_{i,j})\)

\(B=(b_{i,j})\)

\(C=(c_{i,j})\)

\(A\times{B}=(d_{i,j})\)

\(B\times{C}=(e_{i,j})\)

\((A\times{B})\times{C}=(f_{i,j})\)

$A\times{(B\times{C})}=(g_{i,j}) $

由矩阵的乘法得

\(d_{i,j}=a_{i,1}\times{b_{1,j}}+a_{i,2}\times{b_{2,j}}+...+a_{i,n}\times{b_{n,j}}\)

\(e_{i,j}=b_{i,1}\times{c_{1,j}}+b_{i,2}\times{c_{2,j}}+...+b_{i,n}\times{c_{n,j}}\)

\(f_{i,j}=d_{i,1}\times{c_{1,j}}+d_{i,2}\times{c_{2,j}}+...+d_{i,n}\times{c_{n,j}}\)

\(g_{i,j}=a_{i,1}\times{e_{1,j}}+a_{i,2}\times{e_{2,j}}+...+a_{i,n}\times{e_{n,j}}\)

故对任意\(i,j=1,2,...,n\)有,

$f_{i,j}=d_{i,1}\times{c_{1,j}}+d_{i,2}\times{c_{2,j}}+...+d_{i,n}\times{c_{n,j}} $

\(=(a_{i,1}\times{b_{1,1}}+a_{i,2}\times{b_{2,1}}+...+a_{i,n}\times{b_{n,1}})\times{c_{1,j}}+(a_{i,1}\times{b_{1,1}}+a_{i,2}\times{b_{2,1}}+...+a_{i,n}\times{b_{n,1}})\times{c_{2,j}}\)

$+...+(a_{i,1}\times{b_{1,n}}+a_{i,2}\times{b_{2,n}+...+a_{i,n}}\times{b_{n,n}})\times{c_{n,j}} $

\(=a_{i,1}\times{(b_{1,1}\times{c_{1,j}}+b_{1,2}\times{c_{2,j}}+...+b_{1,n}\times{c_{n,j}})}+a_{i,2}\times{(b_{2,1}\times{c_{1,j}}+b_{2,2}\times{c_{2,j}}+...+b_{2,n}\times{c_{n,j}})}\)

$+...+a_{i,n}\times{(b_{n,1}\times{c_{1,j}}+b_{n,2}\times{c_{2,j}}+...+b_{n,n}\times{c_{n,j}})} $

$=a_{i,1}\times{e_{1,j}}+a_{i,2}\times{e_{2,j}}+...+a_{i,n}\times{e_{n,j}}=g_{i,j} $

故\((A\times{B})\times{C}=A\times{(B\times{C})}.\)

证明来自互联网：https://www.cnblogs.com/Jakson/articles/4557558.html

那这又有什么用呢？让我们步入正题。

矩阵快速幂

对于一个\(n×n\)的矩阵，我们称其为方阵，方阵可以进行幂运算，定义为\(C=A^n\)。显然，对于任意一个矩阵，只有它是方阵才能进行幂运算(不然n和m都不相等啊)。

又因为矩阵乘法满足结合律，所以我们可以用快速幂来进行矩阵的幂运算。

矩阵乘法的好处在于它能将有用的状态储存在一个矩阵中，并通过一次乘法得到一次DP的值(DP状态转移方程必须是线性的，虽然有些题目中状态转移方程不是线性的，但是我们可以通过一定的转化将其转化为线性方程)

一般的形式有行向量左乘矩阵和列向量右乘矩阵，例如

\[\begin{bmatrix}f_{i-1}&f_{i}\end{bmatrix}\times\begin{bmatrix}0&1\\1&1\end{bmatrix}=\begin{bmatrix}f_{i}&f_{i+1}\end{bmatrix} \]

实现时一般在向量中补0变为方阵

为什么必须线性？请看下文分解：

我们知道线性的状态转移方程如斐波那契数列

\(f_i=f_{i-1}+f_{i-2}\)

我们能通过单次递推得到斐波那契数列的第n项

它又可以表示为：

\(f_i=1\times{f_{i-1}}+1\times{f_{i-2}}\)

而在进行矩阵乘法时，我们恰恰是用上次的答案来转化出新的状态，所以我们可以用矩阵乘法优化线性DP

结合快速幂：

matrix operator^(int k){
	matrix ans,x=*this;
	memset(ans.g,0,sizeof(ans.g));
	for(int i=1;i<=N;i++)ans.g[i][i]=1;
	for(;k;k>>=1,x=x*x)if(k&1)ans=ans*x;
   return ans;
}

matrix power(matrix a,int b){
   matrix ans;
	memset(ans.g,0,sizeof(ans.g));
	for(int i=1;i<=N;i++)ans.g[i][i]=1;
	while(b){
		if(b&1)ans=multiply(ans,a);
		a=multiply(a,a);
		b>>=1;
   }
  	return ans;
}

关于单位矩阵：单位矩阵相当于乘法中的1，在每次快速幂之前需初始化ans矩阵，代码如下：

memset(ans.g,0,sizeof(ans.g));
for(int i=1;i<=N;i++)ans.g[i][i]=1;

因为使用快速幂，所以矩阵快速幂的时间复杂度基本为\(O(N^3\times{log(n)})\)，其中N为矩阵大小

下面讲解一些例题(讲解可能没有题解中清晰，所以emm，希望更清晰的讲解？点击查看题解)：

  1. 斐波那契数列

矩阵快速幂模板题。

因为

\(f_{i}=1\times{f_{i-1}}+1\times{f_{i-2}}\)

\(f_{i-1}=1\times{f_{i-1}}+0\times{f_{i-2}}\)

所以得到转移矩阵：

\[\begin{bmatrix}1&1\\1&0\end{bmatrix} \]

转移为：

\[\begin{bmatrix}1&1\\1&0\end{bmatrix}\times\begin{bmatrix}f_{i-1}\\f_{i-2}\end{bmatrix}=\begin{bmatrix}f_i\\f_{i-1}\end{bmatrix} \]

用列向量

\[\begin{bmatrix}f_2=1\\f_1=1\end{bmatrix} \]

乘转移矩阵\(n\)次得到\(G_{1,1}=f_{n+1},G_{2,1}=f_{n}\)。输出\(G_{2,1}\)即可。

代码如下：

#include<bits/stdc++.h>
#define int long long
#define inf 0x3f3f3f3f
using namespace std;
int read(){
	int w=0,h=1;char ch=getchar();
	while(ch<\'0\'||ch>\'9\'){if(ch==\'-\')h=-h;ch=getchar();}
	while(ch>=\'0\'&&ch<=\'9\'){w=w*10+ch-\'0\';ch=getchar();}
	return w*h;
}
const int mod=1e9+7;
int n,q;
struct matrix{
	int g[105][105];
	matrix operator*(const matrix&b)const{
		matrix c;
		for(int i=1;i<=2;i++)
			for(int j=1;j<=2;j++){
				c.g[i][j]=0;
				for(int k=1;k<=2;k++){
					c.g[i][j]=(c.g[i][j]+g[i][k]*b.g[k][j]%mod)%mod;
				}
			}
		return c;
	}
	matrix operator^(int k){
		matrix ans,x=*this;memset(ans.g,0,sizeof(ans.g));
		for(int i=1;i<=2;i++)ans.g[i][i]=1;
		for(;k;k>>=1,x=x*x)if(k&1)ans=ans*x;
		return ans;
	}
}f;
signed main(){
	q=read();
	f.g[1][1]=1;f.g[1][2]=1;
	f.g[2][1]=1;f.g[2][2]=0;
	f=f^q;
	printf("%lld\n",f.g[2][1]%mod);
	return 0;
}
  1. \([NOI2012]随机数生成器\)

矩阵快速幂进阶的运用。题目中给出了状态转移方程，根据转移方程列出转移矩阵：

\[\begin{bmatrix}a&0\\c&1\end{bmatrix} \]

初始化\(res\)矩阵:\(res.G_{1,1}=X_0,res.G_{2,2}=1\)；由于我们要推的是\(X_n\)，第一项为\(X_0\)，所以不应该初始化为单位矩阵

发现模数m有点大，用龟速乘防止爆\(longlong\)(什么，你不知道龟速乘？)

#include<bits/stdc++.h>
#pragma GCC optimize(2)
#define int long long
using namespace std
int read(){
    int w=0,h=1;char ch=getchar();
    while(ch<\'0\'||ch>\'9\'){if(ch==\'-\')h=-h;ch=getchar();}
    while(ch>=\'0\'&&ch<=\'9\'){w=w*10+ch-\'0\';ch=getchar();}
    return w*h;
}
struct node{
	int g[5][5];
}f,res,ans;
int m,a,c,x0,n,mod;
void init(）{
	res.g[1][1]=x0%m;res.g[1][2]=1;
	f.g[1][1]=a%m;f.g[1][2]=0;
	f.g[2][1]=c%m;f.g[2][2]=1；
}
int cnt(int x,int y){
	int s=0;
	while(y){
		if(y&1)s=(s+x)%m;
		x=x*2%m;
		y>>=1;
	}
	return s%m;
}
node multiple(node x,node y){
    memset(ans.g,0,sizeof(ans.g))；
    for(int i=1;i<=2;i++)
    for(int j=1;j<=2;j++)
    for(int k=1;k<=2;k++)
    ans.g[i][j]=(ans.g[i][j]+cnt(x.g[i][k]，y.g[k][j]))%m;
    return ans;
}
void Fast(int n）{
	while(n){
		if(n&1)res=multiple(res,f);
		f=multiple(f,f);
		n>>=1;
	}
}
signed main(){
	m=read()；a=read();c=read();x0=read();n=read();mod=read();
	init();
	Fast(n);
	printf("%lld",res.g[1][1]%mod);
    return 0;
}
  1. \(Isaac\)

洛谷爆冷题，谔谔的题面

这是我刚进长郡初中队时教练给我们出的烤试题，然后我当然是不会的，所以就一直放在收藏里，然后某天我点开了这道题，发现居然是矩阵乘法的标签，于是就谔谔开心的做起来，然后我当然还是不会了的，于是就A掉了这题。

看到这题，第一反应当然是图论，但是看到k的要求和数据范围，立刻就懵圈了知道是矩阵乘法，可是图论和矩阵乘法有什么关系呢？\(QwQ\)，于是想了想。

图论和矩阵乘法到底有什么关系呢？请看下文分解：

我们知道：\(C_{i,j}=\sum_{k=1}^{n} A_{i,k}×B_{k,j}\)

而在图论的路径计数中，我们怎样计算从i到j的路径条数呢？与多源最短路算法Floyd相似，我们枚举中转点，根据乘法原理和加法原理，我们得出如下代码：

for(int i=1;i<=N;i++)
	for(int j=1;j<=N;j++)
		for(int k=1;k<=N;k++)
			c[i][j]=c[i][j]+c[i][k]*c[k][j];

与矩阵乘法鲸人的相似！！！

那么，方阵的幂又表示什么呢？

当图G的指数为一时，显然就是一步从i到j的方案数，我们考虑\(G^2\)。

\(G^2\)时，我们枚举了一个中转点，于是原来的路径断成了两部分，而这正是走两步从i到j的方案数！这时，我们就解决了在k时刻恰好到达终点的问题。

但是这题还有一个难点：珂爱的小怪怪。但只要我们抓住这题的关键：数据范围，就能迎刃而解了。

发现珂爱的小怪怪的游走规律不超过4，而1,2,3,4的最小公倍数仅有12，于是我们考虑分情况讨论。我们把每种情况图中怪物将会到达的点标记下来，在计算答案时判断怪物是否在该点并且血量是否足够通过当前边，快速幂前将12种情况乘在一起达到优化效果，注意矩阵乘法没有交换律，所以第0种情况最后乘。由于题目询问最小值，套上二分即可。

对怪物更清晰的讲解见沼泽愕鱼

新鲜出炉的代码：

#include<bits/stdc++.h>
#define int long long
using namespace std;
int read(){
	int w=0,h=1;char ch=getchar();
	while(ch<\'0\'||ch>\'9\'){if(ch==\'-\')h=-1;ch=getchar();}
	while(ch>=\'0\'&&ch<=\'9\'){w=w*10+ch-\'0\';ch=getchar();}
	return w*h;
}
const int mod=10000;
int n,m,St,En,K,NaCly_Fish;
struct node{
	int g[55][55];
	node operator * (node b){
		node c;memset(c.g,0,sizeof(c.g));
		for(int i=1;i<=n;i++){
			for(int j=1;j<=n;j++){
				for(int k=1;k<=n;k++){
					c.g[i][j]|=g[i][k]&b.g[k][j];//二分只需要判断能否到达，所以采用较快的位运算
				}
			}
		}
		return c;
	}
	node operator ^ (int x){
		node ans,t=*this;memset(ans.g,0,sizeof(ans.g));
		for(int i=1;i<=n;i++)ans.g[i][i]=1;
		for(int i=x;i;i>>=1,t=t*t)if(i&1)ans=ans*t;
		return ans;
	}
}a[12],ans,res,check[13],mp;
signed main(){
	n=read();m=read(）;St=read();En=read();K=read();
	for（int i=1;i<=m;i++){
		int u=read(),v=read(),val=read();
		mp.g[u][v]=mp.g[v][u]=val;
	}
	NaCly_Fish=read();//知道为什么是NaCly_Fish吗？因为我做了沼泽愕鱼，沼泽愕鱼里食人鱼数量是NFish……
	for(int i=0;i<12;i++)check[i]=mp;
	for(int i=1;i<=NaCly_Fish;i++){
		int T=read();
		for(int j=0;j<T;j++){
			int cur=read（);
			for(int k=j;k<12;k+=T)
			for(int l=1;l<=n;l++)
			check[k].g[l][cur]=0;//标记在一个循环（轮回）中怪物会到达的点
		}
	}
	int l=0,r=3e9;
	while(l+1<r){
		int mid=(l+r)>>1;
		for(int i=0;i<12;i++)a[i]=check[i]；
		for(int i=0;i<12;i++){
			for(int j=1;j<=n;j++){
				for(int k=1;k<=n;k++){
					if(a[i].g[j][k])a[i].g[j][k]=a[i].g[j][k]<=mid;//判断在一个循环中第i个单位时间到达的点是否有怪且血量是否足够
				}
			}
		}
		memset(res.g,0,sizeof(res.g));
		for(int i=1;i<=n;i++)res.g[i][i]=1;
		for(int i=1;i<12;i++)res=res*a[i];res=res*a[0];
		int t=K/12,step=K%12;//t为12种情况一起乘的次数，step为剩余要乘的情况数
		memset(ans.g,0,sizeof(ans.g));
		for(int i=1;i<=n;i++)ans.g[i][i]=1;
		ans=ans*(res^t)；
		for(int i=1;i<=step;i++)ans=ans*a[i];
		if(ans.g[St][En])r=mid;
		else l=mid;
	}
	if(l>=(1<<31-3))cout<<"\'IMP0SSBLE!!";
	else cout<<r;
	return 0;
}

分类：

技术点：

相关文章：

  • 2021-11-08
  • 2020-06-23
  • 2021-12-03
  • 2021-12-03
  • 2021-12-03
  • 2021-04-04
  • 2021-11-02
  • 2022-01-06
猜你喜欢
  • 2019-08-14
  • 2021-12-03
  • 2021-12-03
  • 2021-12-03
相关资源
相似解决方案
热门标签
Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode