发表更新OI笔记16 分钟读完 (大约2472个字)

矩阵基础

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

一 基础知识

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;
}

我们把里面的a、base、rec假想成一个矩阵，那么rec *= base、base *= 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]; 
}

矩阵基础

https://big-news.cn/2019/07/30/矩阵基础/

作者

ce-amtic

发布于

2019-07-30

更新于

2023-04-15

许可协议

CC BY-NC-SA 4.0

#

评论

Your browser is out-of-date!

Update your browser to view this website correctly.&npsb;Update my browser now

×

×