矩阵基础

又是一个乱堆一气的文章…
索引:矩阵的加减乘运算,矩阵快速幂,矩阵加速线性递推……

一 基础知识

1 概念

一个$n\times m$的数阵被称作矩阵(matrix)。举个栗子,就像下面这样:
$$ \begin{bmatrix} 1&2&3 \newline 4&5&6 \newline 7&8&9 \newline 10&11&12 \end{bmatrix} $$
矩阵可以用大写字母表示,像这样:
$$ A = \begin{bmatrix} 1&2&3 \newline 4&5&6 \newline 7&8&9 \newline 10&11&12 \end{bmatrix} $$
$A_{i,j}$表示矩阵第$i$行第$j$列的元素,如$A_{3,2} = 8$。

主对角线

形如$A_{i,i}$的元素所组成的一条矩阵的对角线被称为矩阵的主对角线。

矩阵的阶

$n$行$m$列的矩阵被称为$n\times m$阶矩阵。如上面的$A$矩阵是$4\times 3$阶矩阵,注意这里的$\times$符号不表示相乘,就是说$3\times 4$和$4\times 3$不是同一个阶数。

2 加减运算

两个阶数相同的矩阵可以进行加减法,就像这样:
若$ A = \begin{bmatrix} 1&2&3 \newline 4&5&6 \end{bmatrix} $,$ B = \begin{bmatrix} 7&8&9 \newline 10&11&12 \end{bmatrix} $
则$A+B = \begin{bmatrix} 1+7 & 2+8 & 3+9 \newline 4+10 & 5+11 & 6+12 \end{bmatrix} = \begin{bmatrix} 8 & 10 & 12 \newline 14 & 16 & 18 \end{bmatrix}$
$A-B = \begin{bmatrix} 1-7 & 2-8 & 3-9 \newline 4-10 & 5-11 & 6-12 \end{bmatrix} = \begin{bmatrix} -6 & -6 & -6 \newline -6 & -6 & -6 \end{bmatrix}$

由上也可知,矩阵的加法是满足交换律的。

3 矩阵相乘

仅有形如这样的两个矩阵可以进行乘法:$n\times k$阶矩阵$A$和$k\times m$阶矩阵$B$。它们的乘积会是一个$n\times m$阶矩阵。

满足相乘条件的矩阵的乘法原则是:若$C = A\times B$,则有$C_{i,j} = \sum\limits_{u = 1}^k A_{i,u}·B_{u,j}$。一个例子:
$$ \begin{bmatrix} 1&2&3 \newline 4&5&6 \end{bmatrix} \times \begin{bmatrix} 1&2 \newline 3&4 \newline 5&6 \end{bmatrix} = \begin{bmatrix} 22 &28 \newline 49&64 \end{bmatrix}$$
具体的运算细节如下图:
矩阵相乘

注意:矩阵乘法是不满足交换律的!


二 矩阵快速幂

1 理论

快速幂能在较短的时间内求出$a^b$的值,是因为巧妙的把$b$进行了二进制分解并划分了求$a^b$的子问题。

那么对于矩阵而言,可不可以用类似于快速幂的办法快速的求出$n\times n$阶矩阵$A$的幂次$A^b$的值呢?
注意:因为矩阵乘法的限制,只有阶数形如$n\times n$的矩阵有幂运算。

显然是可以的,这里有一份快速幂代码:

1
2
3
4
5
6
7
8
9
LL QuickPow(LL a,LL b){ 
LL base = a,rec = 1;
while(b){
if(b & 1) rec *= base;
base *= base;
b >>= 1;
}
return rec;
}

我们把里面的abaserec假想成一个矩阵,那么rec *= basebase *= base我们都可以实现。唯一的问题是rec = 1应该怎么定义?

抛开矩阵,”1”应该满足的条件是:对于任意非零实数$a$,有$a\times 1 = a$。
回到矩阵上来,”1”矩阵$base$应该满足的条件是:对于任意非零矩阵$A$,有$A\times base = A$。

假如$A$为$n\times m$阶,那么显然$base$应该是$m\times m$阶,否则得到的积矩阵不会是$n\times m$阶的。进一步研究发现:$base$应该是一个主对角线为$1$,其余元素都为$0$的矩阵,像这样:$\begin{bmatrix} 1&0&0 \newline 0&1&0 \newline 0&0&1 \end{bmatrix}$(当$m = 3$时)。

2 代码

于是我们得到矩阵快速幂的代码:
模板

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
#define LL long long

const int CN = 101;
const int R = 1e9+7; //结果对 R 取模

class matrix{ //定义一个矩阵
public: LL a[CN][CN]; int n,m; //n*m阶
//matrix() {memset(a,0,sizeof(a)); n=m=0;}
//matrix(int nn,int mm) {memset(a,0,sizeof(a));n = nn; m = mm;}
void MakeOne(){ //构造 "1" 矩阵
for(int i=1;i<=m;i++) a[i][i] = 1;
}
matrix operator *(const matrix& b)const{ //重定义乘法
// n*m x b.n*b.m => n*b.m
// k = m = b.n
matrix rec = (matrix){{{0}},n,b.m}; //乘法运算的答案
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
for(int k=1;k<=b.m;k++)
(rec.a[i][k] += (a[i][j]*b.a[j][k])%R) %= R;
return rec;
}
void operator *=(const matrix& b) {*this = (*this)*b;} //重定义乘法
};
//矩阵快速幂
matrix MatrixQuickPow(matrix &mr,LL k){
matrix rec = (matrix){{{0}},mr.m,mr.m}; //构造 '1' 矩阵
rec.MakeOne(); //构造 '1' 矩阵
while(k){ //类似于普通的快速幂
if(k & 1) rec *= mr;
mr *= mr;
k >>= 1;
}
return mr = rec;
}

时间复杂度大概是$O(k \log_2 b)$,其中$k$是一次矩阵乘法的复杂度,大概是$n^3$级别($n$是阶数)。


三 矩阵加速线性递推

矩阵快速幂可以被用来优化递推。实际上并不是只能优化线性递推,一个非线性递推的例子:请在本站搜索“摆花”。

1 斐波那契数列问题

众所周知的斐波那契数列(Fibonacci Sequence),它的递推式是这样的:$f_1 = f_2 = 1,f_i = f_{i-1} + f_{i-2}(i>2)$。

若求$f_i$项的值,朴素的想法是$O(i)$递推。有没有什么$\log$级别的算法呢?

我们不妨定义一些矩阵:设矩阵$F(i) = \begin{bmatrix} f_i&f_{i-1} \end{bmatrix}$,有没有可能从$F(i-1)$推出$F(i)$呢?

假设存在矩阵$base$,使$F(i-1)\times base = F(i)$。不难发现$F(i-1)\times base^2 = F(i+1)$,于是有一般规律$F(i)\times base^k= F(i+k)$,即可得$F(2)\times base^{i-2}= F(i)$。

我们知道$F(2) = \begin{bmatrix} 1&1 \end{bmatrix}$,而在式子$F(2)\times base^{i-2}= F(i)$中,后面的$base^{i-2}$是可以利用矩阵快速幂快速求得的。于是我们有了一种想法:假如知道$base$,我可以在$\log$级别的时间复杂度里推出$F(i)$,进而推出$f_i$。

$base$所满足的条件是$F(i-1)\times base = F(i)$,即$\begin{bmatrix} f_{i-1}&f_{i-2} \end{bmatrix} \times base=\begin{bmatrix} f_i&f_{i-1} \end{bmatrix}$。至少$base$应该是一个$2\times 2$阶的(请注意推广到一般形式,如果$F(i)$是$1\times n$阶的呢?),否则两矩阵相乘得不到一个$2\times 2$阶矩阵。

我们知道:$f_{i-1}\times 1+f_{i-2}\times 1 = f_i,f_{i-1}\times 1+f_{i-2}\times 0 =f_{i-1}$。即:
$$F(i-1)\times \begin{bmatrix} 1\newline 1 \end{bmatrix} = f_i $$ $$ F(i-1)\times \begin{bmatrix} 1\newline 0 \end{bmatrix} = f_{i-1} $$
故有:$F(i-1)\times \begin{bmatrix} 1&1\newline 1&0 \end{bmatrix} = F(i)$,即$base = \begin{bmatrix} 1&1\newline 1&0 \end{bmatrix}$。

于是利用式子$F(2)\times base^{i-2}= F(i)$,我们就求出$F(i) = \begin{bmatrix} f_i&f_{i-1} \end{bmatrix}$。再求$f_i$就简单不过了。

注意:当递推式不同时,$base$的值是会变的。

代码
模板

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
#define LL long long

const int CN = 101;
const LL R = 1e9+7;

class matrix{ //定义矩阵
public: LL a[CN][CN]; int n,m;
//matrix() {memset(a,0,sizeof(a)); n=m=0;}
//matrix(int nn,int mm) {memset(a,0,sizeof(a));n = nn; m = mm;}
void MakeOne(){
for(int i=1;i<=m;i++) a[i][i] = 1;
}
matrix operator *(const matrix& b)const{
// n*m x b.n*b.m => n*b.m
// k = m = b.n
matrix rec = (matrix){{{0}},n,b.m};
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
for(int k=1;k<=b.m;k++)
(rec.a[i][k] += (a[i][j]*b.a[j][k])%R) %= R;
return rec;
}
void operator *=(const matrix& b) {*this = (*this)*b;}
}F,base;
matrix MatrixQuickPow(matrix mr,LL k){ //矩阵快速幂
matrix rec = (matrix){{{0}},mr.m,mr.m}; //构造 '1' 矩阵
rec.MakeOne(); //构造 '1' 矩阵
while(k){
if(k & 1) rec *= mr;
mr *= mr; k >>= 1;
}
return rec;
}
//主求解函数
LL f(LL i){
if(i <= 2) return 1;
//构造初始矩阵
F = (matrix){{{0}},1,2}; //F(2)
F.a[1][1] = 1; F.a[1][2] = 1; //[f2 f1]
//构造base
base = (matrix){{{0}},2,2};
base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;

F *= MatrixQuickPow(base,i-2);
return F.a[1][1]; //答案是矩阵第一行(一共只有一行...)第一个元素
}

一个很有趣的结论:$ \text{gcd}(f_a,f_b) = f_{\text{gcd}(a,b)} $
例题:LG P1306

2 推广

假如现在不是斐波那契数列,而要求这么一个递推式在第$i$项的值:$f_1=f_2=f_3=1,f_x=f_{x-3}+f_{x-1} (x>3)$,该怎么用矩阵加速呢?

设矩阵$F(i) = \begin{bmatrix} f_{i}&f_{i-1}&f_{i-2} \end{bmatrix}$,则有$F(3) = \begin{bmatrix} f_1&f_2&f_3 \end{bmatrix}$。设$F(i)\times base^k = F(i+k)$,则有$F(i) = F(3)\times base^{i-3}$。通过上面的方法推导出$base = \begin{bmatrix} 1&1&0 \newline 0&0&1 \newline 1&0&0 \end{bmatrix}$(请自己手推一遍,按照上面的分别构造$f_i,f_{i-1},f_{i-2}$的方法!)。

然后求出第$i$项的值就是一个矩阵快速幂+矩阵乘法的事了。

代码
模板

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#define LL long long

const int CN = 101;
const LL R = 1e9+7;

class matrix{ //定义矩阵
/*同上省略*/
}F,base;
matrix MatrixQuickPow(matrix mr,LL k){ //矩阵快速幂
/*同上省略*/
}
//主求解函数
LL f(int i){
if(i <= 3) return 1;//防止快速幂出锅(负数)
//构造初始矩阵
F = (matrix){{{0}},1,3};
F.a[1][1] = F.a[1][2] = F.a[1][3] = 1;
//构造base
base = (matrix){{{0}},3,3};
base.a[1][1] = base.a[1][2] = base.a[2][3] = base.a[3][1] = 1;

F *= MatrixQuickPow(base,i-3);
return F.a[1][1];
}

3 再推广

广义斐波那契数列问题:

广义的斐波那契数列是指形如$f_n=p\times f_{n-1}+q\times f_{n-2}$的数列。今给定数列的两系数$p$和$q$,以及数列的最前两项$f_1$和$f_2$,另给出两个整数$n$和$m$,试求数列的第$n$项$f_n$除以$m$的余数。

这个问题其实依然可以矩阵加速。设矩阵$F(i) = \begin{bmatrix} f_i&f_{i-1} \end{bmatrix}$,有:
$$ F(i-1) \times \begin{bmatrix} p\newline q \end{bmatrix} = f_i$$ $$ F(i-1) \times \begin{bmatrix} 1\newline 0 \end{bmatrix} = f_{i-1}$$ 所以就是$ F(i-1) \times \begin{bmatrix} p&1\newline q&0 \end{bmatrix} = F(i)$。

再套矩阵加速就跟板子一样了。

代码
模板

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
#define LL long long

const int CN = 3;

LL R; //对 R 取余

class matrix{ //定义矩阵
/*同上省略*/
}F,base;
matrix MatrixQuickPow(matrix mr,LL k){ //矩阵快速幂
/*同上省略*/
}
//主求解函数
LL f(LL i,LL f1,LL f2,LL p,LL q){ //f_n = p·f_n-1 + q·f_n-2
if(i == 1) return f1; //特判 1
if(i == 2) return f2; //特判 2
//构造初始矩阵
F = (matrix){{{0}},1,2}; //F(2)
F.a[1][1] = f2; F.a[1][2] = f1; //[f2 f1]
//构造base
base = (matrix){{{0}},2,2};
base.a[1][1] = p; base.a[1][2] = 1; base.a[2][1] = q;

F *= MatrixQuickPow(base,i-2);
return F.a[1][1];
}

$$ - - - - \mathcal{End} - - - - $$


评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×