矩阵乘法学习笔记
【前言】
矩阵乘法通常用于加速线性递推式，能够在较为优秀的时间里完成状态的转移。
题目特征通常有：
- 单次转移规则简单但转移次数很多。
- 转移规则满足结合律。
当然还可以扩展到图上应用，之后都会提到。
【前置芝士】
1.【矩阵】
一个 \(n\times m\) 的矩阵可以看成是一个 \(n\times m\) 的二维数组。
它同样有着和实数一样的某些运算。
2.【矩阵加法】
只有规模相同的矩阵才能相加减，规则为对应位置相加减，即若 \(C=A+B\)，则有：
减法同理。
3.【矩阵乘法】
较为复杂，但同样可以理解。
假设 \(A\) 是 \(n\times m\) 的矩阵，\(B\) 是 \(m\times p\) 的矩阵，\(C=A\times B\)，则有：
显然，只有第一个矩阵的列数等于第二个矩阵的行数时，两个矩阵才可以相乘。
值得思考的是为什么矩阵乘法的定义这么复杂，换句话说，为什么不指定更简单的运算规则达到同样的效果呢？
其实，从矩阵的本质出发，我们很容易找到答案。
矩阵的本质是简化的线性方程组，每一个位置代表了对应未知数的系数。
可以尝试随便写几个线性方程组，简化为矩阵后进行相乘，发现结果确实是上述规则。（具体可以看这里）
同时根据这一点，我们也可以证明矩阵乘法是否符合交换律、结合律、分配率。
结果是：矩阵乘法满足结合律和分配率，不满足交换律。
4.【单位矩阵】
类似于普通乘法中的 \(1\)，设单位矩阵为 \(S\)，有 \(A\times S=A\)。
其构造方式为：
即主对角线上的数字为 \(1\)，其余为 \(0\)。
不难验证其正确性。
【主要思想】
1.【从斐波那契说起】
著名的斐波那契数列有 \(F_i=F_{i-1}+F_{i-2}\) 的递推式，求第 \(n\) 项的时间复杂度为 \(O(n)\)。
但是如果结合矩阵乘法，会有意想不到的优化。
2.【构造状态矩阵】
要得到 \(F_i\)，显然只需要知道 \(F_{i-1},F_{i-2}\) 即可，所以每次只需要保留最后两位。
设 \(F(n)\) 表示一个 \(1\times 2\) 的矩阵，\(F(n)=[Fib_n\ Fib_{n+1}]\)。
我们希望根据 \(F(n-1)=[Fib_{n-1}\ Fib_{n}]\) 得到 \(F(n)\)，这需要另一个矩阵的辅助。
3.【构造转移矩阵】
不难发现我们需要将 \(F(n-1)\) 的第一项变为 \(Fib_{n-1}+Fib_{n}\)，第二项变为 \(Fid_{n}\)。
所以不难构造出一个转移矩阵：
根据上述乘法原则，不难发现 \(F(n)=F(n-1)\times A\)。
4.【矩阵快速幂加速递推】
假设我们需要第 \(F(n)\)，设 \(F(0)=[0\ 1]\)，显然 \(F(n)=F(0)\times A^n\)。
所有具有乘法结合律的乘方运算都可以用快速幂来加速，矩阵同理。
不仅如此，快速幂是贯穿矩阵学习的重要思想，请不了解的同学寻找相关资料学习。
5.【Code】
代码大概长酱紫：
void Mul(int f[2],int a[2][2]){
int c[2];
memset(c, 0, sizeof(c));
for(int j=0;j<2;j++)
for(int k=0;k<2;k++)
c[j] = (c[j] + 1LL * f[k] * a[k][j] % MOD) % MOD;
memcpy(f, c, sizeof(c));
}
void Mulself(int a[2][2]){
int c[2][2];
memset(c, 0, sizeof(c));
for(int i=0; i<2; i++)
for(int j=0; j<2; j++)
for(int k=0; k<2; k++)
c[i][j] = (c[i][j] + 1LL * a[i][k] * a[k][j] % MOD) % MOD;
memcpy(a, c, sizeof(c));
}
int f[2]={0,1};
int a[2][2]={{0,1},{1,1}};
for(; n; n >>= 1){
if(n & 1) Mul(f, a);
Mulself(a);
}
printf("%lld\n", f[0]);
6.【小结】
通过这道题，我们对矩阵乘法加速递推有了初步了解。
一般来说，此类题目有以下特点：
- 可以抽象出一个长度为 \(n\) 的一维向量，该向量在每个单位时间内只发生一次变化。
- 变化形式为线性递推。（加上若干数或乘上一个数）
- 该递推式规则本身保持不变。（根据这个构造转移矩阵）
- 递推轮数很多但向量长度 \(n\) 不大。
那么我们可以考虑矩阵乘法优化，这个长度为 \(n\) 的矩阵称为状态矩阵，上面出现过的 \(A\) 称为转移矩阵。
解题的关键就是构造出合适的转移矩阵和状态矩阵。
转移矩阵一般的构造规则是：
如果状态矩阵中的第 \(x\) 个数对下一时间单位的第 \(y\) 个数有影响，那么将 \(A_{x,y}\) 赋为恰当的值。
时间复杂度为 \(O(n^3\log T)\)，\(T\) 为递推总轮数。
【经典例题】
【典例一】
已知一个数列 \(a\)，它满足：
\[a_x= \begin{cases} 1 & x \in\{1,2,3\}\\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases} \]求 \(a\) 数列的第 \(n\) 项。
不难想到推导过程。
可以定义状态矩阵：\(F(n)=[a_n,a_{n-1},a_{n-2}]\)。
显然 \(F(n+1)=[a_n+a_{n-2},a_n,a_{n-1}]\)。
所以有：
【典例二】
给定大小为 \(n\) 的正整数可重集 \(S\)，每次操作可以加入数 \(a+b(a,b\in S)\)，求 \(k\) 次操作后 \(S\) 的和的最大值。
显然每次贪心的寻找最大的两个数相加，关键是求和。
这需要状态矩阵加上一维表示和。
令 \(F(n)=[f(n-1),f(n),S(n)=\sum_{i=1}^n f(i)]\)（\(f(n),f(n-1)\) 分别表示当前集合内的最大值和次大值）
显然 \(F(n+1)=[f(n),f(n-1)+f(n),S(n)+f(n-1)+f(n)]\)。
则有：
【典例三】
给定矩阵 \(A\) 和 \(k\)，求 \(A+A^2+A^3+...+A^k\)。
关于矩阵加法的定义已经给出，不再赘述。
直接计算的复杂度为 \(O(n^3k)\)，显然不可取。
考虑分治，每次将序列二分并分治下去，计 \(F(k)=A+A^2+A^3+...+A^k\)。
显然对于每个 \(k\)：
将矩阵乘法在结构体中重载，代码还挺好写的，之后都采用这种方式。
struct Matrix{
int a[50][50];
int MOD = 10;
Matrix() {memset(a, 0, sizeof(a));}
Matrix operator * (const Matrix &b){
Matrix c;
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
for(int k=0; k<n; k++)
c.a[i][j] = (c.a[i][j] + a[i][k] * b.a[k][j] % MOD) % MOD;
return c;
}
Matrix operator + (const Matrix &b){
Matrix c;
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
c.a[i][j] = (a[i][j] + b.a[i][j]) % MOD;
return c;
}
Matrix operator ^ (int k){
Matrix c, b;
for(int i=0; i<n; i++) c.a[i][i]=1;//这里用到了单位矩阵
memcpy(b.a, a, sizeof(a));
for(; k; k >>= 1){
if(k & 1) c = c * b;
b = b * b;
}
return c;
}
} A;
Matrix dfs(int k){
if(k == 1) return A;
Matrix P = dfs(k >> 1);
P = P + P * (A ^ (k >> 1));
if(k & 1) P = P + (A ^ k);
return P;
}
【典例四】
太长了，自己看
同样将问题转化为一维向量。
定义长度为 \(n\times m+1\) 的状态矩阵 \(F\)，下标为 \(0\) 到 \(n\times m\)。
记 \(num(i,j)=(i-1)\times m+j\)。
用 \(F[num(i,j)]\) 记录格子 \((i,j)\) 中的石头数量。
特别的，令 \(F[0]=1\)，把它当做得到石头的窗口，原因后文介绍。
设 \(F(k)\) 为 \(k\) 秒后的状态矩阵，显然 \(F(0)=[1,0,0,...,0]\)。
可惜本题的转移矩阵是不断变化的，而且每个格子还不一样。
但是可以发现 \(lcm(1,2,3,4,5,6)=60\)，所以我们可以记录下前 \(60\) 个状态矩阵。
设 \(A_k\) 表示第 \(k\) 个，构造方法为：
- 若第 \(k\) 秒 \((i,j)\) 的操作字符为 \(N\)，则令 \(A_k[num(i,j),num(i-1,j)]=1\)。其余同理。
- 若第 \(k\) 秒 \((i,j)\) 的操作字符为数字 \(x\)，则令 \(A_k[0,num(i,j)]=x\)。
- 令 \(A_k[0,0]=1\)。
- 令其余为 \(0\)。
注意 \(F[0]=1\)，显然以上矩阵可以胜任我们的状态转移工作。
设 \(A=\Pi_{i=1}^{60} A_i,t=60\times k+b\)，那么 \(F(t)=F(0)\times A^k\times \Pi_{i=1}^b A_i\)，那么：
【典例五】
询问有向图上从 \(s\) 点恰好经过 \(k\) 个点到达 \(t\) 点的方案数
经典的图上问题转化为矩阵问题。
首先将图转化为邻接矩阵（这注定了 \(n\) 的范围不会很大，这也是此类题目的标志之一）
令 \(C=A\times A\)，那么 \(C[i][j]=\sum A[i][k]×A[k][j]\) ，实际上就等同于 \(i\) 到 \(j\) 恰好经过两条边的路径数量。
\(k\) 步同理。
其实相当于 Floyd 算法的矩阵乘法实现。
而利用快速幂相当于 Floyd 算法的倍增实现。
【典例六】
给定一张 \(T\) 条边的无向连通图，求从 \(S\) 到 \(E\) 经过 \(N\) 条边的最短路长度。
从路径个数到最短路长度，其实没变多少。
矩阵乘法时将乘法变为取 min，其实这更像 Floyd 了。
主要思想依然是用矩阵来优化。
【典例七】
太长了，自己看
看到 \(n\leq 100\) 和这个转移方式，矩阵不难想到。
显然根据异或的性质，异或偶数次等于没有异或。
所以我们只需要统计每个点 \(u\) 对 \(1\) 号点的影响次数，如果影响次数为偶数次就异或上，然后得到答案。
我们可以用一个 \(0/1\) 矩阵每个位置 \((i,j)\) 表示节点 \(i\) 是否对节点 \(j\) 有贡献，利用矩阵快速幂加速计算。
可惜这样我们每次都要乘一次，时间复杂度：\(O(qn^3\log a)\)，好像不太行。
考虑倍增优化。
假设初始矩阵为 \(A\)，我们可以预处理出 \(A,A^2,A^4,..A^{2^{32}}\) 相当于对其进行了二进制拆分。
时间复杂度变为 \(O(n^3\log a+qn^2\log a)\)，好像可以。
【总结】
根据以上例题的学习总结，我们完成了矩阵从基本运算到数列递推到图上问题的完美转化。
值得注意的是，根据典例二、四，我们发现在矩阵乘法加速递推遇到求和或常数项时，通常有一下方法：
在状态矩阵中添加一个额外位置，始终储存数列和或常数项 \(1\)，并乘上转移矩阵中的适当系数。
通过这种方法将系数累加到状态矩阵的相应位置。
根据典例六我们发现，利用矩阵乘法求解 Floyd 时并不需要像 Floyd 那样先枚举中间节点。
因为这本质上是矩阵乘法，而普通的 Floyd 算法基于 DP 实现，有严格的阶段性。
至此矩阵乘法入门的全部内容已整理完毕。
引用资料&特别鸣谢：
- 矩阵乘法优化DP复习。
- 【学习笔记】动态规划—矩阵递推加速。
- 《算法竞赛进阶指南》。
完结撒花
