矩阵加速学习笔记
_Jocularly_
·
·
算法·理论
矩阵
矩阵运算
矩阵其实类似于一个二维数组，我们通常定义一个矩阵用大写字母来表示，比如 A = \begin{bmatrix}2 & 1 \\ 1 & 2 \\ \end{bmatrix}，那么我们在此基础上，可以对矩阵进行计算，有的时候会起到意想不到的功能。
矩阵加减
加减比较好理解，两个相同大小的矩阵加减就是对应位置加减，例如 A = \begin{bmatrix}2 & 1 \\ 1 & 2 \\ \end{bmatrix}，B = \begin{bmatrix}1 & 2 \\ 2 & 1 \\ \end{bmatrix}，A - B = \begin{bmatrix}2 - 1 & 1 - 2 \\ 1 - 2 & 2 - 1 \\ \end{bmatrix} = \begin{bmatrix}1 & -1 \\ -1 & 1 \\ \end{bmatrix}，相加是相同的。
矩阵数乘
顾名思义，就是矩阵乘上一个常数，法则为每个位置对应乘上这个数，例如 A = \begin{bmatrix}2 & 1 \\ 1 & 2 \\ \end{bmatrix}，A \times 2 = \begin{bmatrix}2\times 2 & 1\times 2 \\ 1\times 2 & 2\times 2 \\ \end{bmatrix} = \begin{bmatrix}4 & 2 \\ 2 & 4 \\ \end{bmatrix}。
矩阵互乘
矩阵互乘，又称矩阵乘法，是我们利用来解决问题的基本方式。首先给出一个标准定义的计算法则：
两个大小分别为 $m \times n$ 和 $n \times p$ 的矩阵 $A, B$ **相乘**的结果为一个大小为 $m \times p$ 的矩阵。将结果矩阵记作 $C$，则：$C_{i,j} = \sum_{k = 1}^{n} A_{i,k} \times B_{k,j}$。
看上去有点匪夷所思，但是如果找找规律就能发现，我们可以类比二维数组来看，定义二元组 $(i,j)$ 表示矩阵 $A$ 的第 $i$ 行和矩阵 $B$ 的第 $j$ 列，这一行和列对应位置相乘求和，放在数组下标为 $(i,j)$ 的位置，得到的数组就为矩阵乘法的结果。
举个例子，还是这两个矩阵 $A = \begin{bmatrix}2 & 1 \\ 3 & 4 \\ \end{bmatrix}$，$B = \begin{bmatrix}5 & 6 \\ 8 & 7 \\ \end{bmatrix}$，那么定义 $A \times B = C$，$C = \begin{bmatrix}2\times5 + 1\times 8 & 2\times 6 + 2\times 7 \\ 3\times 5 + 4\times 8 & 3\times 6 + 4\times 7 \\ \end{bmatrix}$。
这里就要注意，我们要求二元组 $(i,j)$ 对应相乘，也就是说长度要是相等的，也就是要求第一个矩阵的长一定要等于第二个矩阵的宽。
### 矩阵乘法计算法则
#### 交换律
矩阵乘法不满足交换率，一定要记住这一点，因为 $A \times B$ 和 $B \times A$ 的结果并不相同，**不可以进行交换**。
#### 结合律
包括两部分，满足加括号或删括号。加括号指 $ABC = (AB)C = A(BC)$，删括号指 $A(B + C) = AB + AC$。
### 矩阵幂运算
一个大小为 $n \times n$ 的矩阵 $A$ 可以与自身进行乘法，得到的仍是大小为 $n \times n$ 的矩阵，记作 $A^2 = A \times A$。进一步地，还可以递归地定义任意高次方 $A^k = A \times A^{k - 1}$，或称 $A^k = \underbrace{A \times A \times \cdots \times A}_{k \text{ 次}}$。
### 单位矩阵
就和单位 $1$ 一样，在普通乘法中，乘上一个单位 $1$ 结果是不变的。类似的，这个乘上任意矩阵都等于原矩阵的矩阵称之为单位矩阵，定义单位矩阵为 $I = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}$。
## 矩阵快速幂
矩阵快速幂就是矩阵乘法加上快速幂，因为我们可以使用类似于快速幂的思想快速求出矩阵幂运算，详细内容看代码。
```cpp
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
#define mod 1000000007
int n,k;
struct matrix{//矩阵
int data[105][105];
matrix(){
memset(data,0,sizeof(data));//初始化
}
void init(){//单位矩阵
for(int i=1;i<=n;i++){
data[i][i] = 1;
}
}
};
matrix a;
matrix ans;
matrix mul(matrix a,matrix b){//矩阵乘法核心
matrix res;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
for(int k=1;k<=n;k++){
res.data[i][j] += (a.data[i][k] * b.data[k][j]) % mod;
res.data[i][j] %= mod;
}
}
}
return res;
}
matrix quick_power(int k){//快速幂
matrix cnt;
cnt.init();
while(k){
if(k & 1) cnt = mul(cnt,a);
a = mul(a,a);
k >>= 1;
}
return cnt;
}
signed main(){
cin >> n >> k;
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
cin >> a.data[i][j];
}
}
ans = quick_power(k);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
cout << ans.data[i][j] << " ";
}
cout << endl;
}
return 0;
}
```
## 矩阵加速
### 例题
看一道[例题](https://www.luogu.com.cn/problem/P1939)：
ax=
\begin{cases}
1 & x \in{1,2,3}\
a{x-1}+a_{x-3} & x \geq 4
\end{cases}
求 $a$ 数列的第 $n$ 项对 $10^9+7$ 取余的值。
看起来和矩阵乘法驴唇不对马嘴，但是我们可以把它转化为矩阵乘法。我们要求的是 $a_x$，为了求出下面的几项，我们还需要求出 $a_{x-1}$ 和 $a_{x-2}$，那么我们把它抽象成矩阵，得到 $\begin{bmatrix} a_{x} \\ a_{x-1} \\ a_{x-2}\end{bmatrix}$，采用题目中的递推式子，推出可以求出矩阵的式子，把矩阵化为 $\begin{bmatrix} a_{x} = a_{x-1} \times1 + a_{x-2} \times 0 + a_{x-3} \times 1 \\ a_{x-1} = a_{x-1} \times1 + a_{x-2} \times 0 + a_{x-3} \times 0 \\ a_{x-2} = a_{x-1} \times0 + a_{x-2} \times 1 + a_{x-3} \times 0\end{bmatrix}$，那么我们根据矩阵乘法的逆操作，把他转化为 $\begin{bmatrix} a_{x} \\ a_{x-1} \\ a_{x-2}\end{bmatrix} = \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0\end{bmatrix} \times \begin{bmatrix} a_{x-1} \\ a_{x-2} \\ a_{x-3}\end{bmatrix}$，那么提取出来一个 $x-1$，就可以转化为矩阵快速幂来求解。
定义初始矩阵 $A = \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0\end{bmatrix}$，那么 $\begin{bmatrix} a_{x} \\ a_{x-1} \\ a_{x-2}\end{bmatrix} = A^{x-1} \times \begin{bmatrix} 1 \\ 1 \\ 1\end{bmatrix} = A^{x-1}$。此时我们求得答案即为第一行第一个。详见代码
```cpp
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
#define mod 1000000007
int t;
int n;
struct matrix{
int data[5][5];
matrix(){
memset(data,0,sizeof(data));
}
void init(){
for(int i=1;i<=3;i++){
data[i][i] = 1;
}
}
void base(){
data[1][1] = 1;
data[1][3] = 1;
data[2][1] = 1;
data[3][2] = 1;
}
};
matrix ans;
matrix mul(matrix a,matrix b){
matrix res;
for(int i=1;i<=3;i++){
for(int j=1;j<=3;j++){
for(int k=1;k<=3;k++){
res.data[i][j] += (a.data[i][k] * b.data[k][j]) % mod;
res.data[i][j] %= mod;
}
}
}
return res;
}
matrix quick_power(int k){
matrix cnt;
matrix a;
cnt.init();
a.base();
while(k){
if(k & 1) cnt = mul(cnt,a);
a = mul(a,a);
k >>= 1;
}
return cnt;
}
signed main(){
cin >> t;
while(t--){
cin >> n;
ans = quick_power(n-1);
cout << ans.data[1][1] << endl;
}
return 0;
}
```
------
那么会了这道题，我们再看一道[例题](https://www.luogu.com.cn/problem/P1962)，题意非常的清楚，即为求斐波那契数列第 $n$ 项，那么此时我们根据上一题的算法，列出需要求的内容为 $\begin{bmatrix} F_{i} \\ F_{i-1} \end{bmatrix}$，此时根据题目中的递推式得到 $\begin{bmatrix} F_{i} \\ F_{i-1} \end{bmatrix} = \begin{bmatrix} 1\times F_{i-1} + 1\times F_{i-2} \\ 1\times F_{i-1} + 0 \times F_{i-2} \end{bmatrix}$，定义转换 $base$ 矩阵为 $A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$，那么我们可以得到公式为 $\begin{bmatrix} F_{i} \\ F_{i-1} \end{bmatrix} = A \times \begin{bmatrix} F_{i-1} \\ F_{i-2} \end{bmatrix} = A^{i-1}$，直接计算即可。
```cpp
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define endl "\n"
#define mod 1000000007
int n;
struct matrix{
int data[5][5];
matrix(){
memset(data,0,sizeof(data));
}
void init(){
for(int i=1;i<=2;i++){
data[i][i] = 1;
}
}
void base(){
data[1][1] = 1;
data[1][2] = 1;
data[2][1] = 1;
}
};
matrix ans;
matrix mul(matrix a,matrix b){
matrix res;
for(int i=1;i<=2;i++){
for(int j=1;j<=2;j++){
for(int k=1;k<=2;k++){
res.data[i][j] += (a.data[i][k] * b.data[k][j]) % mod;
res.data[i][j] %= mod;
}
}
}
return res;
}
matrix quick_power(int k){
matrix cnt;
matrix a;
cnt.init();
a.base();
while(k){
if(k & 1) cnt = mul(cnt,a);
a = mul(a,a);
k >>= 1;
}
return cnt;
}
signed main(){
cin >> n;
ans = quick_power(n-1);
cout << ans.data[1][1];
return 0;
}
```
## 广义矩阵乘法
定义 $\otimes$ 有交换律为 $a \otimes b = b \otimes a$，定义 $\oplus$ 有结合律为 $a \oplus b \oplus c = a\oplus(b\oplus c)$，称 $\otimes$ 对 $\oplus$ 有分配率为 $a\otimes (b \oplus c) = (a\otimes b) \oplus (a\otimes c)$。
我们发现狭义矩阵乘法是 $C_{i,j} = \sum_{k = 1}^{n} A_{i,k} \times B_{k,j}$。此时我们发现在原矩阵乘法中乘法对加法有分配率，那么此时我们可以拓展至广义矩阵乘法，只要 $\otimes$ 满足交换律结合律，且 $\otimes$ 对 $\oplus$ 有分配率，那么 $C_{i,j} = \bigoplus_{k = 1}^{n} A_{i,k} \otimes B_{k,j}$，且该矩阵运算满足结合律，但不满足交换律。发现 $\max$ 或 $\min$ 对常见运算均满足以上条件，所以广义矩阵乘法通常运用在优化动态规划。
例如最大子段和，状态转移方程为 $dp_{i} = \max(dp_{i-1} + a_{i},a_{i})$，那么此时我们发现 $\max$ 和 $+$ 满足广义矩阵乘法，那么我们可以将转移方程转化为矩阵乘法，$\begin{bmatrix} dp_i \\ 0 \end{bmatrix} = \begin{bmatrix} \max(dp_{i-1}+a_{i},a_i+0) \\ 0 \end{bmatrix} = \begin{bmatrix} dp_{i-1} \\ 0 \end{bmatrix} \times \begin{bmatrix} a_i & a_i \\ -\infty & 0 \end{bmatrix}$，这里的 $-\infty$ 相当于狭义乘法中的 $0$，相当于不可能取这个，当为 $\min$ 时取 $+\infty$。