矩阵类优化 dp——矩阵加速优化

Idea

了解矩阵优化之前，我们不妨先来说说什么是矩阵。

如其名，矩阵就是按照长方形阵列排列的实数（也可以是复数）集合，类似于这样：

$$S=\begin{bmatrix}1&2&3\\4&5&8\end{bmatrix}$$

而我们在优化 dp 转移时主要用到的就是矩阵的乘法。

设我们有两个矩阵 $S,T$，若 $S$ 和 $T$ 之间可以作乘法，首先需要保证矩阵 $S$ 的列数和矩阵 $T$ 的行数相同。所以不妨设 $S$ 的规模为 $x\times y$，$T$ 的规模为 $y\times z$，若有 $S\times T=R$，则 $R$ 的规模为 $x\times z$，且对于 $R$ 有：

$$R_{i,j}=\sum_{k=1}^yS_{i,k}\times T_{k,j}$$

那么问题来了，为什么我们要把好端端的一个递推式变成一个矩阵进行转移呢？答案很简单，因为矩阵乘法可以用矩阵快速幂快速求解——矩阵乘法是满足结合律的！

也就是说，对于三个矩阵 $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_i=f_{i-1}+f_{i-2}$$

也就是说，我们的转移需要用到前面两项的值，我们也只关心这两项。于是我们可以将初始矩阵状态设为最前面的两项，即 $[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_1,f_2]\times\begin{bmatrix}1&1\\0&1\end{bmatrix}=[f_2,f_3]$$

扩展一下，若我们想得到向量 $[f_n,f_{n+1}]$，则可以有递推式：

$$[f_1,f_2]\times\begin{bmatrix}1&1\\0&1\end{bmatrix}^{n-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$ 也是需要维护的。一共十一个变量。

所以我们的初始矩阵就为：

$$[a_1,b_1,c_1,a_2,b_2,c_2,1,1,1,w,z]$$

转移一步之后的矩阵为：

$$[a_2,b_2,c_2,a_3,b_3,c_3,1,2,4,w^2,z^2]$$

可以推出转移矩阵为：

$$\begin{bmatrix} 0&0&0&q&0&0&0&0&0&0&0\\ 0&0&0&0&v&0&0&0&0&0&0\\ 0&0&0&0&0&y&0&0&0&0&0\\ 1&0&0&p&1&1&0&0&0&0&0\\ 0&1&0&1&u&1&0&0&0&0&0\\ 0&0&1&1&1&x&0&0&0&0&0\\ 0&0&0&1&0&2&1&1&1&0&0\\ 0&0&0&t&0&1&0&1&2&0&0\\ 0&0&0&r&0&0&0&0&1&0&0\\ 0&0&0&0&1&0&0&0&0&w&0\\ 0&0&0&0&0&1&0&0&0&0&z\\ \end{bmatrix}$$

简单来说就是，将变量设进向量里，将常量放进转移矩阵里，推出前两步后的初始矩阵向量后，通过他们之间的递推关系推出转移矩阵即可。

注意到本题模数最大为 $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_i=\sum_{j=\max(0,i-m)}^{i-1}f_j$$

开始时将 $f_0$ 的初值设为 $1$ 即可，表示开头一定会放上一块木板。

而观察到 $m$ 很小，转移阶段又很大（$10^{500}$ 过于魔鬼），所以可以考虑矩阵加速优化转移。

但是显然还不够，因为 $10^{500}$ 次方显然过于庞大了，甚至还需要用到高精，且 $g$ 的拆分数目足足可以达到 $2^{|s|}$ 种，这是不好的。考虑到 $g$ 数组的拆分方式十分奇怪，思考一下题目中是否有什么厉害的性质可以让我们优化。

注意我们转移时用的是矩阵，而矩阵所满足的一个性质是：

$$A^{n+m}=A^m\times A^n$$

正确性是很显然的。所以例如 $g(123)$，我们就可以将它拆解为：

$$f(1+2+3)+f(12+3)+f(1+23)+f(123)=f(1+2)\times f(3)+f(12)\times f(3)+f(1)\times f(23)+f(123)$$

考虑递推求解 $g$，设 $g'_i$ 表示 $s$ 前 $i$ 位数字组成的数字串的答案（我们最终要求出 $g'_n$），则显然当我们算到 $g(123)$ 时，$g(1)$ 和 $g(12)$ 都已经算过了。于是可以将上面的式子进一步化简为：

$$(f(1+2)+f(12))\times f(3)+f(1)\times f(23)+f(123)=f(123)+g'_1\times f(23)+g'_2\times f(3)$$

所以扩展地来看，只要我们知道每一个后缀数字串的 $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 提高组] 魔法值

甚至状态和转移方程题目中都直接给出来了：

$$f_{u,i}=\text{xor}_{v}f_{v,i-1}$$

将 $u$ 之外的节点与 $u$ 的连边状态显式的写出来，设 $e_{u,v}$ 表示 $u$ 与 $v$ 之间是否连边，则有：

$$f_{u,i}=\text{xor}\ f_{v,i-1}e_{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$ 为合法决策，则有：

$$dp_i=\max(dp_j+c_i)$$

而 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$ 次方矩阵，按照十进制拼接起来即可。

P4569 [BJWC2011] 禁忌

AT2371 [AGC013E] Placing Squares

---

Ending

总结下来，矩阵优化 dp 的特征其实非常鲜明，只要看完数据范围基本就可以确定此题是否可以矩阵优化。

具体来说，当我们观察数据范围时，发现：

1. 如果 dp 式子可以写成矩阵，且矩阵规模较小。
2. 转移阶段过大，正常复杂度无法承担。

那么我们就完全可以考虑矩阵优化 dp。

而矩阵优化 dp 的主要用途无非两种：

1. 优化转移递推式。
2. 优化边权较小图上路径统计问题。

列矩阵时，基本步骤如下：

1. 详细列出题目中的 dp 式子，判断广义矩阵乘法的形式。

2. 找出 dp 式子中的变量与常量，并溯源每个变量的求解递推式，将变量列进初始矩阵，常量列进转移矩阵，矩阵加速幂优化转移即可。

以及转移过程中的一些 trick：

1. 若题目中有一些关键节点会改变转移矩阵，考虑到向量乘矩阵和矩阵乘矩阵的不同复杂度，可以考虑先预处理出转移矩阵的二进制整数幂，再对于每个关键节点“暂停转移”，依次处理。

2. 若发现模数超过 long long 范围，记得使用龟速乘。

3. 如果转移矩阵都为 $0/1$ 构成，则可以考虑使用 bitset 优化。

4. 若图上的转移涉及多个时间戳，则可以类似于拆点将每一个时间戳的该点都列出转移矩阵。

5. 阶段较大时，十进制转二进制会非常麻烦，则可以考虑预处理十进制下 $0-9$ 的转移矩阵，在十进制下直接转移即可。