矩阵类优化 dp——矩阵加速优化
Idea
了解矩阵优化之前，我们不妨先来说说什么是矩阵。
如其名，矩阵就是按照长方形阵列排列的实数（也可以是复数）集合，类似于这样：
而我们在优化 dp 转移时主要用到的就是矩阵的乘法。
设我们有两个矩阵 \(S,T\)，若 \(S\) 和 \(T\) 之间可以作乘法，首先需要保证矩阵 \(S\) 的列数和矩阵 \(T\) 的行数相同。所以不妨设 \(S\) 的规模为 \(x\times y\)，\(T\) 的规模为 \(y\times z\)，若有 \(S\times T=R\)，则 \(R\) 的规模为 \(x\times z\)，且对于 \(R\) 有：
那么问题来了，为什么我们要把好端端的一个递推式变成一个矩阵进行转移呢？答案很简单，因为矩阵乘法可以用矩阵快速幂快速求解——矩阵乘法是满足结合律的！
也就是说，对于三个矩阵 \(S,T,R\)，是有 \((S\times T)\times R=S\times(T\times R)\) 的。
这个其实很好证明，可以试着自己证明一下。（基本原理其实是加法满足对于乘法的分配律，且乘法满足结合律）
所以当转移递推式固定时，我们可以通过类似于快速幂的方式，设置好初始矩阵和转移矩阵后，通过矩阵快速幂将转移时关于阶段的一维复杂度优化为 log 级别。具体设置方式可以看下面例题。
值得扩展的一点时，注意只要我们的矩阵运算满足结合律就可以进行快速幂，而乘法和加法相结合的矩阵运算之所以会满足结合律，其实就是因为乘法本身具有结合律且加法对于乘法具有分配律。所以广义地来讲，只要我们的矩阵运算法则中，内层的运算具有结合律，外层的运算对于内层的运算具有分配律，且状态转移递推式恒定，再结合数据范围发现递推方程涉及到底前缀状态极少，转移阶段又十分大时，就可以考虑通过矩阵快速幂进行优化转移。（例如 min 和 +，or 和 and）
设矩阵规模为 \(n\times n\)，转移阶段为 \(m\)，则一般时间复杂度为 \(O(n^3\log m)\)。
Extension
1. 图上 dp 优化
考虑有经典问题：
给定一个图，边权均为 \(0\) 或 \(1\)，求从起点到终点的最短路径。
首先显然可以考虑通过 \(\text{dijkstra}\) 来正常做，复杂度应为 \(O(m\text{ log }m)\).
但由于边权的特殊性，我们完全可以建立一个双端队列，遇到 \(1\) 加队尾，\(0\) 加队首，跑一遍 \(\text{bfs}\)，即可得到答案，从而去掉一个 \(\text{log}\)，时间复杂度 \(O(m)\).
扩展成一个计数问题：
给定一个图，对于两个点之间要么不连边，要么边权为 \(1\)，求从起点出发长度为 \(k\) 到达终点的路径方案数（\(k\) 极大，\(n\le 500\)）。
考虑如果没有括号里的限制，显然可以用动态规划来做。
设 \(dp_{i,j}\) 表示路径长度为 \(i\)，当前达到结点 \(j\) 的路径数目。则初始设 \(dp_{0,start}\) 为 \(1\)，因为边权的特殊性，通过相邻边进行转移即可。时间复杂度 \(O(mk)\)，空间复杂度 \(O(m)\).
而 \(k\) 是极大的，\(n\) 的范围却是允许立方级别复杂度，自然就想到了矩阵加速优化。通过相临边的关系构造转移矩阵，用矩阵快速幂转移即可。时间复杂度 \(O(n^3 \text{ log } k)\)，空间复杂度 \(O(n^2)\)。
若起点终点不固定，则改变初始化方式和最终统计方式即可，算法流程不变。
Problem
1. P1962 斐波那契数列
一道很基础的矩阵加速幂优化，就拿这道题来谈谈实现细节。
首先我们有一个固定的递推式，即：
也就是说，我们的转移需要用到前面两项的值，我们也只关心这两项。于是我们可以将初始矩阵状态设为最前面的两项，即 \([f_1,f_2]\).
首先来考虑这个矩阵怎么变成 \([f_2,f_3]\)。设 \([f_1,f_2]\times S=[f_2,f_3]\)，则我们都知道 \(S\) 的哪些信息呢？
首先 \(S\) 的规模一定是 \(2\times 2\) 的，这个毋庸置疑。故设 \(S=\begin{bmatrix}a&b\\c&d\end{bmatrix}\)，考虑转移。
由矩阵乘法的规则可知，\(f_2=f_1\times a+f_2\times c\)，所以 \(a=0,c=1\)；又由于 \(f_3=f_1\times b+f_2\times d\)，所以 \(b=1,d=1\)。
即我们有：
扩展一下，若我们想得到向量 \([f_n,f_{n+1}]\)，则可以有递推式：
此时输出向量的第一个元素即可。
矩阵加速幂函数如下：
void ksm(int k){//矩阵a的k次方
while (k){
if (k&1) ans=ans*a;
a=a*a;
k>>=1;
}
return ans;
}
重载矩阵运算符代码如下：
struct node{
ll s[N][N];//每一个结构体都代表一个矩阵
node(){
memset(s,0,sizeof(s));//自动初始化
}
}a,ans;
node operator *(const node&a,const node&b){//重载运算符板子
node zc;
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
for (int k=1;k<=3;k++)
zc.s[i][j]=(zc.s[i][j]+a.s[i][k]*b.s[k][j]%mod)%mod;
return zc;
}
一种更好理解的方式是，我们可以将转移矩阵 \(S\) 中第 \(i\) 行第 \(j\) 列的元素 \(S_{i,j}\) 看作是要将向量的第 \(i\) 个元素乘上 \(S_{i,j}\) 后赋值给新向量的第 \(j\) 个元素。注意矩阵的运算位置不能乱改。
2. P1707 刷题比赛
这个递推式和数据范围一看就很矩阵加速。
不过递推式涉及到的东西有点多，让我们来好好捋一捋。
首先 \(a,b,c\) 在转移过程中会互相涉及，所以我们显然是要同时转移这些变量的。
那么当我们的转移阶段为 \(k+2\) 时，都有哪些东西需要知道呢？
首先 \(a_{k+1},b_{k+1},c_{k+1},a_k,b_k,c_k\) 肯定都是需要知道的，常数 \(1\) 也需要用到，还需要一个自增的 \(k\)（每次转移时都从 \(1\) 那里给他一个 \(1\)），以及一个自增的 \(k^2\)，\(w^k\) 和 \(z^k\) 也是需要维护的。一共十一个变量。
所以我们的初始矩阵就为：
转移一步之后的矩阵为：
可以推出转移矩阵为：
简单来说就是，将变量设进向量里，将常量放进转移矩阵里，推出前两步后的初始矩阵向量后，通过他们之间的递推关系推出转移矩阵即可。
注意到本题模数最大为 \(1e16\)，直接做乘法应该会爆 long long，所以需要使用龟速乘。
时间复杂度 \(O(11^3\log^2 n)\).
3. P4159 [SCOI2009] 迷路
给定一个有向图，边权为 \(1\) 到 \(9\)，问从起点到终点长度为 \(k\) 的路径数（\(k\le 10^9\)，\(n\le 10\)）。
注意到 \(n\) 极小，\(k\) 极大，考虑转化为上面的经典问题进行解决。
发现和上面经典问题唯一的区别就是边权不再只有 \(1\) 了，而是变成了 \(1\) 到 \(9\)，但依然很小。于是我们很容易就想到，将边拆开，拆成边权 \(-1\) 个点，对于每条边从起点到终点依次串联起拆成的点，如此以来，可以使得新图中边权均为 \(1\)，就转化为了上述经典问题。
可惜如此以来，\(n\) 相当于最大变为了 \(1000\) 左右，由于时间复杂度为 \(O(n^3\text{ log } k)\)，显然是过不了大数据的。
考虑优化。
首先肯定是要拆的，拆边肯定是不行的，而图上只有点和边，所以就考虑拆点。
考虑对于每一个点，向周围连的边边权均为 \(1\) 到 \(9\)，所以我们可以将每一个点拆为 \(9\) 个点，分别表示周围边权为 \(1\) 到周围边权为 \(9\)。对于拆的点与外界的点，若与外界连的边边权为 \(i\)，则从代表边权为 \(i\) 的点与其代表边权为 \(1\) 的点连边；对于拆点内部，代表边权为 \(i\) 的点向代表边权为 \(i+1\) 的点连边。如此，可以保证当前图与原图是等价的。于是愉快地跑一边矩阵加速即可！
由于拆的节点最多只有 \(90\) 个，显然可以通过此题。
由此我们可以发现，拆点其实就是将原图改成一个与之等价的图，使得题目中的特殊限制被满足，从而达到可以套板子解决原问题的目的。而对于拆点的方式，可以考虑将边拆为点，也可以考虑将原点拆为若干个点，但应在满足题意的情况下，尽量拆更少的点，使得新图更加简单，从而降低复杂度。
4. P2151 [SDOI2009]HH去散步
给定一个无向图，边权均为 \(1\)，求起点到终点长度为 \(k\) 的路径数目（\(n\le 50\)，\(m\le60\)，\(k\le2^{30}\)）。特殊限制是上一步走过的路下一步不能重复走。
经典问题的痕迹很明显，通过我们上一题的经验之谈，要做的就是通过拆点使得新图满足特殊限制。
由于边数其实也不多，就可以考虑将每个点拆成相邻边数个点，第 \(i\) 个点代表从第 \(i\) 条相邻的边走过来的点，然后向该点除了第 \(i\) 条边之外的边所连的点连一条单向边即可构造出新图。
但还没完。考虑对于起始点，可以向周围任意一个点前进，所以我们还可以额外设置一个 \(0\) 号点，向起始点周围的所有点连边，然后将 \(dp_{0}\) 的初始值设为 \(1\) 即可跑矩阵加速。
对于每一条边，都会割出两个点，最多割出 \(120\) 个点，时间复杂度 \(O(m^3\text{ log } k)\)，实测效率较为优秀（单题最优解前五）。
由这道题其实也可以看出，割点其实就是将每个点的不同状态放大化，使其“人格分裂”，从而将【节点的转移】转化为【状态的转移】，从而达到实现满足特殊条件约束的目的。
5. P3176 [HAOI2015]数字串拆分
首先我们发现 \(f\) 数组是可以通过 dp 求出的。
不妨将拆分看作是给 \(n\) 个物品插板，两个木板之间的物品个数就是拆分出来的一个数。设 \(f_i\) 为在第 \(i\) 个物品后面插板之后所能得到的拆分数目，则显然，上一个木板可以放在往前的 \(m\) 个空隙当中的任何一个。即有转移方程：
开始时将 \(f_0\) 的初值设为 \(1\) 即可，表示开头一定会放上一块木板。
而观察到 \(m\) 很小，转移阶段又很大（\(10^{500}\) 过于魔鬼），所以可以考虑矩阵加速优化转移。
但是显然还不够，因为 \(10^{500}\) 次方显然过于庞大了，甚至还需要用到高精，且 \(g\) 的拆分数目足足可以达到 \(2^{|s|}\) 种，这是不好的。考虑到 \(g\) 数组的拆分方式十分奇怪，思考一下题目中是否有什么厉害的性质可以让我们优化。
注意我们转移时用的是矩阵，而矩阵所满足的一个性质是：
正确性是很显然的。所以例如 \(g(123)\)，我们就可以将它拆解为：
考虑递推求解 \(g\)，设 \(g'_i\) 表示 \(s\) 前 \(i\) 位数字组成的数字串的答案（我们最终要求出 \(g'_n\)），则显然当我们算到 \(g(123)\) 时，\(g(1)\) 和 \(g(12)\) 都已经算过了。于是可以将上面的式子进一步化简为：
所以扩展地来看，只要我们知道每一个后缀数字串的 \(f\) 值，就可以转移 \(g'\)，得到 \(g\) 的答案了！
而为了避免高精，设 \(f'_{i,j}\) 表示 \(s\) 第 \(i\) 位数字到第 \(j\) 位数字组成的 \(f\) 的转移矩阵，则通过 \(f'_{i,j}=f(s_i)^{10^{n-i}}\times f'_{i+1,j}=f(10^{n-i})^{s_i}\times f'_{i+1,j}\) 即可推出 \(f_{i,j}'\)，所以预处理出所有 \(f(10^i)\) 即可。其中 \(s_i\) 表示代表 \(s\) 中第 \(i\) 个数字的矩阵。
最后输出 \(g'_n\) 即可。
6. CF576D Flights for Regular Customers
题目关心的只有走过的边数，相当于每条边的边权都为 \(1\)。很熟悉的感觉。
又观察到数据范围。更熟悉了。矩阵加速的味道十分浓重。
如果将当前走过的边数看作时间戳，可以先思考一下若每条边的时间戳都为 \(0\)，即每条边一开始就存在时，该如何做。
显然，跑一边 bfs 即可解决。无解的情况很好判断，只要 \(1\) 号节点在边全都存在的情况下不能到达 \(n\) 号节点，则在任意时间戳上皆无法达到。
那么有解的时候呢？我们设 \(dp_{i,j}\) 表示走到节点 \(j\) 时，长度为 \(i\) 的路径是否存在，则显然 \(dp_{i,j}\) 可以通过 \(j\) 号节点的所有相邻节点 \(k\) 的 \(dp_{i-1,k}\) 或运算过来。由于递推式子恒定，所以可以考虑矩阵乘法，并且用 与和或 代替 乘和加。
那有 \(d_i\) 该怎么做呢？考虑到将所有的 \(d_i\) 排序，则两个 \(d_i\) 之间的转移方程显然是不变的，我们在每到一个时间戳的时候停止转移，然后从目前所有能到的地方跑一边向终点跑一边多源 bfs，将 bfs 到的最短距离与目前已经走过的边数相加即为当前的答案，修改转移矩阵之后再继续转移即可。注意由于每个时间戳边数不同，所以先到达终点时的答案不一定最优，需要将每个时间戳上的答案全部算出来，取最小值就是答案。
时间复杂度 \(O(n^3\log |d|\times m)\)，好像过不了。
观察到 \(dp\) 值只有 \(0\) 和 \(1\)，所以用 bitset 优化转移即可，时间复杂度会乘上 \(\frac{1}{w}=\frac{1}{32}\)，差不多能过。
注意用 bitset 的写法需要将两个矩阵对齐，大概可以这么写：
node operator *(const node&x,const node&y){
node z;
for (int k=1;k<=n;k++)
for (int i=1;i<=n;i++)
if (x.a[i][k]) z.a[i]|=y.a[k];
return z;
}
7. P6569 [NOI Online #3 提高组] 魔法值
甚至状态和转移方程题目中都直接给出来了：
将 \(u\) 之外的节点与 \(u\) 的连边状态显式的写出来，设 \(e_{u,v}\) 表示 \(u\) 与 \(v\) 之间是否连边，则有：
所以我们可以设置一个 xor 和 乘法 的广义矩阵乘法。但是这个广义矩阵乘法是否满足条件了？
前面提到，广义矩阵乘法的设置需要满足两个条件：
-
内层运算要满足结合律——乘法显然满足 √
-
外层运算要对于内层运算满足分配律——异或对于乘法满足吗？
答案显然是否定的，反例随处可见，直觉上也并不符合。
那怎么办？？这可不好转化啊！
别着急，别忘了我们的 \(e\) 只有 \(0/1\) 两种取值！在 \(0/1\) 状态下的矩阵乘法是否满足呢？\(0\) 就消失，\(1\) 就保留，当然满足。
此时时间复杂度为 \(O(n^3 \log |a|q)\)，还是不太行。
注意到一个细节，向量乘矩阵的时间复杂度其实可以是 \(O(n^2)\) 的，而矩阵乘矩阵则需要 \(O(n^3)\)。我们在求解矩阵快速幂的时候就是因为要边算向量乘矩阵、边算矩阵乘矩阵，所以矩阵乘的时间复杂度才是 \(O(n^3)\)。如果我们将两者分开呢？
观察到我们是可以快速计算一个矩阵的 \(2\) 的整数次幂的，如果我们将所有的整数次幂矩阵都预处理出来，那么在询问的时候，我们就可以通过经典套路二进制拆分，只通过向量乘矩阵得到最终的答案向量。
故最终时间复杂度 \(O(n^3\log|a|+qn^2|a|)\)。
启发性地讲，对于类似于 凑 \(k\) 的多次询问或需要多次用到 凑 \(k\)，且 求解二进制整数幂答案的时间复杂度较低，则可以考虑倍增+二进制拆分 求解。
8. P5059 中国象棋
注意到每一行都答案其实是相对独立且相同的，所以只需要算出一行的答案，即可得到整个棋盘的答案。
我们都关注哪些信息呢？阶段显然是设置成当前处理到了前 \(i\) 个格点，除此之外，我们显然还关心当前格点有没有放置卒且当前行一共有几个卒，这样才能使我们转移出的方案满足题目中给定的两个条件。
那么现在其实就有两种设置状态的选择：
- 设 \(dp_i\) 表示处理完前 \(i\) 个格点的方案数，且第 \(i\) 个格点强制放置时的方案数，此时我们只需要枚举上一次放置的格点，注意一些细节即可同时满足两个条件。
- 设 \(dp_{i,0/1,0/1/2}\) 表示处理完了前 \(i\) 个格点，且第 \(i\) 个格点 没有放置/放置棋子，本行目前 没有棋子/有一个棋子/棋子个数达到要求 时的方案数。则转移时考虑上一个格子的状态即可。
两个 dp 状态设置的区别在于一个强制选择而另一个不强制选择，正常来说应该是第一个 dp 更方便一些的，但是观察到此题的数据范围，\(N\le 10^{18}\)，果断考虑矩阵加速优化，那么由于矩阵转移的特殊性，转移时决策显然是越少越好的，第一种的决策涉及到了前面的所有状态，显然无法写成矩阵，于是我们便采取了第二种状态设计方式。
对于每一个阶段，一共有六种不同的状态，所以我们设计一个 \(6\times 6\) 的矩阵即可优化转移。最后将答案 \(s\) 作一个 \(n\) 次幂即可。
注意到模数的取值范围，记得使用龟速乘。
时间复杂度 \(O(6^3\log N\log P)\).
9. P3597 [POI2015] WYC
看到特殊边权，按照套路，先考虑拆点，但是我们显然不能将起点或者终点为拆得的点的路径算进答案，于是考虑换种思路。
由于 \(3\) 实在是太小了，所以不妨将其直接放入矩阵，设 \(dp_{i,j}\) 表示的是到达节点 \(j\) 的长度为 \(i\) 的路径的条数，而 \(dp_i\) 是可以由 \(dp_{i-1},dp_{i-2},dp_{i-3}\) 转移过来的，于是我们将三维都设进矩阵，根据边权转移矩阵即可。除此之外，为了快速的计算小于等于当前矩阵阶段的路径个数，我们可以再开一维矩阵表示当前的总和，每次转移时将新产生的路径条数全部加到这一维度当中即可。
最后计算时考虑到矩阵快速幂倍增的特殊形式，可以考虑用倍增代替二分，算出最终的答案。
不过注意此题十分毒瘤，个别数据最后可能会爆 long long，所以建议使用 int128 解决本题。
10. P6772 [NOI2020] 美食家
有了前面几道题的铺垫，这道题显得就不是那么难了。先考虑没有美食节该怎么做。
首先依然采用上一题的套路，按照边权写成一个 \(5n\times 5n\) 规模的矩阵。若设 \(dp_i\) 表示在当前阶段到达 \(i\) 的最大预约值，\(j\) 为合法决策，则有：
而 max 满足结合律，加法具有对于 max 的分配律，所以可以由此定义广义矩阵乘法。
那么加上美食节该怎么做呢？
考虑对于每次美食节中断一次转移，特殊转移一遍美食节矩阵之后，再继续正常转移，但这样每两次美食节之间的转移都是 \(O(n^3\log T)\) 的，总时间复杂度就是 \(O(n^3k \log T)\)。
所以考虑使用 T7 的套路预处理出每一个二进制整数次幂的矩阵，两次美食节之间转移时直接调用即可，由于向量乘矩阵的优秀复杂度，总时间复杂度就被我们优化到了 \(O(n^3 \log T+n^2k\log T)\)。
但是注意转换矩阵乘法规则之后不要思维定势，矩阵赋值项设为 \(1\) 而非 \(0\)，无用项设为不可能值（-INF）而非 \(0\)，矩阵的初始化设为不可能值（-INF）而非 \(0\)，以及一些细节问题。
Exercise
P2886 [USACO07NOV]Cow Relays G 将最短路和矩阵加速相结合，用 min + 替换 + *
P1397 [NOI2013] 矩阵游戏 行内转移，行间转移，预处理出一位的 \(0-9\) 次方矩阵，按照十进制拼接起来即可。
AT2371 [AGC013E] Placing Squares
Ending
总结下来，矩阵优化 dp 的特征其实非常鲜明，只要看完数据范围基本就可以确定此题是否可以矩阵优化。
具体来说，当我们观察数据范围时，发现：
- 如果 dp 式子可以写成矩阵，且矩阵规模较小。
- 转移阶段过大，正常复杂度无法承担。
那么我们就完全可以考虑矩阵优化 dp。
而矩阵优化 dp 的主要用途无非两种：
- 优化转移递推式。
- 优化边权较小图上路径统计问题。
列矩阵时，基本步骤如下：
-
详细列出题目中的 dp 式子，判断广义矩阵乘法的形式。
-
找出 dp 式子中的变量与常量，并溯源每个变量的求解递推式，将变量列进初始矩阵，常量列进转移矩阵，矩阵加速幂优化转移即可。
以及转移过程中的一些 trick：
-
若题目中有一些关键节点会改变转移矩阵，考虑到向量乘矩阵和矩阵乘矩阵的不同复杂度，可以考虑先预处理出转移矩阵的二进制整数幂，再对于每个关键节点“暂停转移”，依次处理。
-
若发现模数超过 long long 范围，记得使用龟速乘。
-
如果转移矩阵都为 \(0/1\) 构成，则可以考虑使用 bitset 优化。
-
若图上的转移涉及多个时间戳，则可以类似于拆点将每一个时间戳的该点都列出转移矩阵。
-
阶段较大时，十进制转二进制会非常麻烦，则可以考虑预处理十进制下 \(0-9\) 的转移矩阵，在十进制下直接转移即可。
