学习笔记 - 矩阵加速

似乎是我没有听讲……矩阵加速一直没有学会 QwQ
寒假将至,还是学一学


『预备』

矩阵加速主要是用到了矩阵快速幂,那么显然我们要知道什么是矩阵乘法。
(Tab: $a\times b$ 表示 $a$ 行 $b$ 列)
首先两个矩阵要满足“前一个的列数和后一个的行数相同”才能相乘,不妨设矩阵 $A$ 为 $r\times l$ 的矩阵, $B$ 为 $l\times c$ 的矩阵。那么它们的乘积 $C=AB$ 为一个 $r\times c$ 的矩阵,且:

举个例子:

那么我们可以进入下一步了。


『矩阵快速幂』

由于矩阵乘法要求“前一个的列数和后一个的行数相同”,那么这就意味着在大多数情况下不存在交换律,不过所幸它满足结合律~

这就意味着我们能够像做一般的乘法一样对矩阵做快速幂,那么求 $A^n$ (此时 $A$ 的行数等于列数等于 $x$)的复杂度为 $O(x^3log_2n)$ ,至于那个 $x^3$ 是矩阵乘法的时间复杂度,一般看成常数。

那么我们只需要定义一个矩阵结构体(MATRIX),再重定义它的乘法运算,套用快速幂的版就可以完成矩阵快速幂了。

1
2
3
4
5
6
7
8
9
MATRIX QPow(MATRIX A,int n){
MATRIX ret=one; //(!)
while(n){
if(n&1) ret=ret*A;
A=A*A;
n>>=1;
}
return ret;
}

但是还有一个问题……大家应该注意到了,上面代码中有一个”one”。在快速幂的模板中,储存答案的变量最初要赋值为 “1”,那么对于矩阵,“1”就是单位矩阵——一个从左下到右上的对角线上全是 $1$,其余全是 $0$ 的一个正方形矩阵(这种矩阵乘上任何一个其他矩阵后仍然得到那一个矩阵)。


『矩阵加速』

矩阵加速是用来优化递推的(至少我现在学到的是这样)。

对于递推式 $f(n)=k_1f(n-1)+k_2f(n-2)+…+k_mf(n-m)$ ($k_1,k_2,…,k_m$ 均为常数),当任何一个 $n$ 都满足这个式子(不存在特殊情况),我们就可以通过矩阵加速求得 $f(n)$。

结合几道例题来讲可能会方便一些~

「例题Ⅰ- HihoCoder 1143」

本质上就是求斐波那契数列的第n项,稍微打下表就可以发现。这道题的 $n\leq 10^8$ 按道理来说 $1s$ 可以直接递推……但是既然我们在学矩阵加速,我们还是用矩阵加速的方法。

矩阵加速最重要的就是找到“矩阵乘法的转移式”;($f(n)$ 即第 $n$ 项)首先我们知道 $f(n)=f(n-1)+f(n-2)$,那么求斐波那契第 $n$ 项的转移式就是:

令函数 $g(n)=\begin{smallmatrix}\begin{bmatrix}f(n)\\f(n-1)\end{bmatrix}\end{smallmatrix}$,常量 $M=\begin{smallmatrix}\begin{bmatrix}1 1\\1 0\end{bmatrix}\end{smallmatrix}$ ,那么递推式就可以写成 $g(n)=M\times g(n-1)$ 。

现在再来看 $M$ 和 $g(n)$ 是如何求出来的:
a) 首先我们必须要能够从 $g(n-1)$ 推导到 $g(n)$ 。根据 ①式,我们采取的方法是“提取出等号左边与 $n$ 相关的值并将它排成一列作为 $g(x-1)$”,即将 $f(n-1),f(n-2)$ 排成一列,作为 $g(n-1)$,那么我们就可以得到 $g(n)$。
b) 假设 $g(n)$ 的行数为 $R$,$M$ 一般来说是 $R\times R$ 的一个矩阵。$M$ 的第 $i$ 行代表的是 $g(n)$ 的第 $i$ 行如何由 $g(n-1)$ 转移过来。

这下面就是我推导 $M$ 的方法:

写出 $R$ 个等式,第 $i$ 个等式的左边是 $g(n)$ 的第 $i$ 行;把 $g(n-1)$ 横向写在每个等式的右边,得到

这样就能够看出 $M_{(1,1)}=1;M_{(1,2)}=1;M_{(2,1)}=1;M_{(2,2)}=0$ ,对应的就是矩阵 $M$。

「例题Ⅱ- HihoCoder 1151」

先思考出状压DP的方程式(这里就不推导了)。

$dp[S][i]$ 表示将前 $i-1$ 列填满时第 $i$ 列被填充的状态为 $S$ (具体见HihoCoder的提示)
则有:
$dp[000][i]=dp[111][i-1]$
$dp[001][i]=dp[110][i-1]$
$dp[010][i]=dp[101][i-1]$
$dp[011][i]=dp[111][i-1]+dp[100][i-1]$
$dp[100][i]=dp[011][i-1]$
$dp[101][i]=dp[010][i-1]$
$dp[110][i]=dp[111][i-1]+dp[001][i-1]$
$dp[111][i]=dp[000][i-1]+dp[011][i-1]+dp[110][i-1]$

根据之前的思路,构造出

可能需要多理解一会,由于要将构造过程完整地写出来会占用大量篇幅,就不写了,如果有疑问记得 email ~


『源代码』

两道例题……以供参考。(把MATIRX打包了一个结构体,可以单独剪出来用)

「例题Ⅰ- HihoCoder 1143」

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
/*Lucky_Glass*/ 
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int R=3,C=3,MOD=19999997;
struct MATRIX{
int r,c;
int mat[R][C];
MATRIX operator *(const MATRIX ano)const{
MATRIX ret,A,B;
if(c==ano.r) A=*this,B=ano;
else if(r==ano.c) A=ano,B=*this;
else exit(0);
for(int i=0;i<A.r;i++)
for(int j=0;j<B.c;j++){
ret.mat[i][j]=0;
for(int k=0;k<A.c;k++)
ret.mat[i][j]=int((1ll*ret.mat[i][j]+1ll*A.mat[i][k]*B.mat[k][j])%(1ll*MOD));
}
ret.r=A.r;ret.c=B.c;
return ret;
}
MATRIX QPow(int k){
MATRIX A=*this,ret;
ret.r=ret.c=2;
ret.mat[0][0]=0;ret.mat[0][1]=1;
ret.mat[1][0]=1;ret.mat[1][1]=1;
while(k){
if(k&1) ret=ret*A;
A=A*A;
k>>=1;
}
return ret;
}
};
int n;
int main(){
scanf("%d",&n);
MATRIX one;
one.r=one.c=2;
one.mat[0][0]=0;one.mat[0][1]=1;
one.mat[1][0]=1;one.mat[1][1]=1;
MATRIX fir;
fir.r=1;fir.c=2;
fir.mat[0][0]=0;fir.mat[0][1]=1;
fir=fir*(one.QPow(n));
printf("%d\n",fir.mat[0][0]);
return 0;
}

「例题Ⅱ- HihoCoder 1151」

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
/*Lucky_Glass*/
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll MOD=12357ll,N=100000000ll;
ll n;
/*
0 1 3 7 5 2 6 4
i 000 001 011 111 101 010 110 100
ma- xo xx xx xx xx xo xo xo
tr- xo xo xx xx xo xx xx xo
ix xo xo xo xx xx xo xx xx
111 110 111 011 010 101 001 011
100 110 111
000

0 1 2 3 4 5 6 7
0 | 0 0 0 0 0 0 0 1 dp[0] dp[7]
1 | 0 0 0 0 0 0 1 0 dp[1] dp[6]
2 | 0 0 0 0 0 1 0 0 dp[2] dp[5]
3 | 0 0 0 0 1 0 0 1 dp[3] dp[4]+dp[7]
4 | 0 0 0 1 0 0 0 0 * dp[4] = dp[3]
5 | 0 0 1 0 0 0 0 0 dp[5] dp[2]
6 | 0 1 0 0 0 0 0 1 dp[6] dp[1]+dp[7]
7 | 1 0 0 1 0 0 1 0 dp[7] dp[0]+dp[3]+dp[6]
*/
const int R=8,C=8;
struct MATRIX{
int r,c;
int mat[R][C];
//A's col = B's rol
MATRIX operator *(const MATRIX ano)const{
MATRIX ret,A,B;
if(c==ano.r) A=*this,B=ano;
else exit(0);
for(int i=0;i<A.r;i++)
for(int j=0;j<B.c;j++){
ret.mat[i][j]=0;
for(int k=0;k<A.c;k++)
ret.mat[i][j]=int((1ll*ret.mat[i][j]+1ll*A.mat[i][k]*B.mat[k][j])%(1ll*MOD));
}
ret.r=A.r;ret.c=B.c;
return ret;
}
MATRIX QPow(int k){
MATRIX A=*this,ret;
ret.r=ret.c=8;
for(int i=0;i<8;i++)
for(int j=0;j<8;j++)
if(i==j) ret.mat[i][j]=1;
else ret.mat[i][j]=0;
while(k){
if(k&1) ret=ret*A;
A=A*A;
k>>=1;
}
return ret;
}
};
MATRIX one,fir;
int main(){
scanf("%lld",&n);
fir.r=one.r=one.c=8;fir.c=1;
fir.mat[7][0]=1;
for(int i=0;i<8;i++) one.mat[i][7-i]=1;
one.mat[3][7]=one.mat[6][7]=one.mat[7][3]=one.mat[7][6]=1;
one=one.QPow(n);
fir=one*fir;
printf("%d\n",fir.mat[7][0]);
return 0;
}

The End

Thanks for reading!

Email: lucky_glass@foxmail.com ,欢迎提问~

0%