首页 / 资讯中心 / 文章详情

矩阵快速幂专题:从线性递推到状态转移的建模与实战

矩阵快速幂专题:从线性递推到状态转移的建模与实战 ★ FEATURED ARTICLE
带这期LCY算法进阶培训班之前我把历年矩阵快速幂专题的课件重新翻了一遍。说实话这个专题在校队训练计划里属于“看着简单讲透很难”的类型矩阵乘法人人都学过快速幂也几乎是基础班标配可把它们凑在一起后很多同学会在转移矩阵上一卡就是一晚上。LCY培训班的本期主题定位很准不贪多只围绕线性递推和状态转移这两类模型把矩阵快速幂拆成建模、构造、调参三部分来讲。本文就把这次专题训练里我觉得最值得记录的思路、模板和坑位都整理出来适合刚学完快速幂、想在竞赛里稳定拿分的同学参考。如果你已经会写矩阵快速幂模板重点看转移矩阵怎么构造这部分才是拿分的关键。1. 为什么矩阵快速幂会成为必学专题1.1 这个专题总被放进进阶班的原因校队训练里矩阵快速幂一般不会在基础班讲。原因不是代码难写而是前置知识凑齐需要时间你至少得先懂快速幂还得习惯矩阵乘法的“行乘列”视角更重要的是得能把一个递推式抽象成状态向量和转移矩阵。基础班很多同学还在用循环做斐波那契当n到1e8就当场裂开这时候讲矩阵快速幂思维跨度太大。但放到进阶班就非常合适。矩阵快速幂的知识点封闭逻辑链短两节课能把模板和典型题讲完训练效果立竿见影。更重要的是它能承接很多后续专题概率DP里的马尔可夫转移、AC自动机上的计数DP、图上长度为k的路径统计这些都要用到同一个套路。带训练这几年我发现凡是这期专题认真跟下来的同学后面学状态压缩DP和生成函数时明显轻松很多。所以LCY培训班把它放在进阶班第一周是一个性价比很高的安排。1.2 矩阵快速幂到底解决了什么问题一句话概括凡是能用“当前的若干个状态量线性组合出下一步状态”的问题都可以尝试矩阵快速幂。最经典的场景是k阶线性递推比如斐波那契。朴素循环算第n项是O(n)n到1e18就完全不可行矩阵快速幂把复杂度压到O(k^3 log n)k是状态维度通常只有3到10log n大概是60这个复杂度在竞赛题里可以轻松过。除递推外还有三类问题也常用它。第一是图上的路径计数邻接矩阵的k次幂直接给出两点之间恰好走k步的方案数这个性质在做计数题时几乎是秒杀。第二是概率期望问题很多随机过程的“下一步概率只取决于当前分布”本质上就是行向量乘以转移矩阵。第三是带有前缀和或多项式附加项的递推把S(n)、n、n²这样的量也塞进状态向量里就能统一处理。这几种场景在省选和区域赛里反复出现所以专题训练从来不缺素材。2. 核心原理与转移矩阵构造2.1 快速幂的二进制思想怎么迁移到矩阵普通快速幂求a^k用到的核心观察是幂次k可以拆成二进制比如求2^100只需要算2^64、2^32、2^4再乘起来而不需要连乘100次。每次把一个数平方指数翻倍相当于把问题规模对数缩减。你手头有1块钱每轮允许翻倍要凑出64块并不需要翻倍64次只需要6轮翻倍这就是快速幂的生活化类比。矩阵快速幂完全一样只是把“一个数”换成了“一个方阵”。关键前提是矩阵乘法满足结合律所以A^(ab) A^a × A^b成立。很多人担心矩阵乘法不满足交换律其实这里用不到交换律快速幂本身只依赖结合律。单位矩阵相当于数字里的1任何矩阵乘单位矩阵都等于它自身所以快速幂的初始结果矩阵设为E而不是全零矩阵。这一点看起来基础但真有人在这里翻车后面排查表里还会提。2.2 斐波那契转移矩阵的手推全过程从最经典的递推说起。F(n) F(n-1) F(n-2)已知F(1)1F(0)0求F(n)在模1e97下的值。构造思路是不要只看单个F(n)把相邻两项打包成一个状态向量v(n) [F(n), F(n-1)]^T我们希望从v(n-1)推到v(n)也就是[F(n), F(n-1)]^T M × [F(n-1), F(n-2)]^T矩阵M的每一行对应新状态的一个分量。第一行F(n) 1×F(n-1) 1×F(n-2)所以第一行是[1,1]第二行F(n-1) 1×F(n-1) 0×F(n-2)所以第二行是[1,0]。于是M [[1,1], [1,0]]写代码时套快速幂答案是(M^(n-1) × v(1))的第一个分量其中v(1) [F(1), F(0)]^T [1,0]^T。很多人算出来的F(n)总差一项问题多半出在幂次多乘了一次或少乘了一次。我第一次带训时就强调一个验证方法先用n5手算M^4再乘v(1)得到的第一个值必须是5。如果第一个值是8说明你多乘了一次如果是3说明少乘了一次。2.3 n²这类多项式项怎么塞进矩阵真正的难点来了。如果递推式里带着n²、n、常数比如a(n) 2×a(n-1) 3×a(n-2) n²你会发现单独一个a(n-1)、a(n-2)不够了n²不能写成a(n-1)的常系数组合。技巧在于把n²也当作状态变量用“下一步的n²可以由上一步的n²、n和1线性表示”来构造矩阵。因为从n-1走到n时n² (n-1)² 2(n-1) 1所以只要把(n-1)²、n-1、1都放进向量里就能精确表达。定义旧状态v(n-1) [a(n-1), a(n-2), (n-1)², n-1, 1]^T新状态v(n) [a(n), a(n-1), n², n, 1]^T。转移矩阵的每一行分别写第一行a(n) 2×a(n-1) 3×a(n-2) (n-1)² 2×(n-1) 1所以是[2,3,1,2,1]第二行a(n-1)直接保留所以是[1,0,0,0,0]第三行n² (n-1)² 2×(n-1) 1所以是[0,0,1,2,1]第四行n (n-1) 1所以是[0,0,0,1,1]第五行1始终是1所以是[0,0,0,0,1]这样拼出来的5×5矩阵配合初始向量一次自乘就能算到很大步数。我自己的经验是遇到带多项式附加项的问题先列“从k-1到k每个分量怎么变”再填矩阵行比硬背模板靠谱得多。常数的本质是“步数无关的固定项”所以单独占一个状态位n的幂次需要几阶就放几个连续幂次进去直到它们能互相线性表示。3. 模板代码和性能优化实录3.1 一份可直接上手的C竞赛模板矩阵快速幂的模板代码很短但细节决定成败。我在训练课上给过一份内部模板核心如下#include bits/stdc.h using namespace std; const long long MOD 1000000007LL; struct Mat { static const int SZ 5; int n; long long a[SZ][SZ]; Mat(int _n 0, bool unit false) : n(_n) { memset(a, 0, sizeof(a)); if (unit) { for (int i 0; i n; i) a[i][i] 1; } } Mat operator*(const Mat other) const { Mat res(n); for (int i 0; i n; i) { for (int k 0; k n; k) { if (a[i][k] 0) continue; for (int j 0; j n; j) { res.a[i][j] (res.a[i][j] a[i][k] * other.a[k][j]) % MOD; } } } return res; } }; Mat qpow(Mat base, long long e) { Mat res(base.n, true); while (e 0) { if (e 1) res res * base; base base * base; e 1; } return res; }两个位置需要特别解释。第一res(base.n, true)传的第二个参数是unit表示初始化为单位矩阵这直接对应“快速幂的结果初始为1”这个习惯。第二乘法循环顺序我写的是i-k-j而不是最直观的i-j-k原因下一节细说。使用这个模板求斐波那契第n项时只需要构造2×2的M和初始向量v然后取qpow(M, n-1) × v的第一维。静态数组的好处是避免vector的动态堆分配在反复乘法的场景下会稳定不少。3.2 矩阵乘法的循环顺序与取模优化矩阵乘法三重循环的顺序看起来只是ijk的排列问题实际会影响常数。我实测下来i-k-j的顺序在绝大多数情况下比i-j-k快尤其在矩阵乘法被快速幂反复调用时。原因是i-k-j的写法把a[i][k]固定在相对外层可以在内层循环里连续使用同时other.a[k][j]的访问更规整缓存命中率也会好一些。我在代码里还加了一句if (a[i][k] 0) continue这个对稀疏矩阵很有效比如邻接矩阵、转移矩阵中有大量0时能省掉不少无意义计算但如果矩阵是稠密的加这个判断反而多花时间要根据题目类型取舍。取模也要注意。MOD是1e97时两个long long相乘大约是1e18在long long范围内不会溢出但累加多轮后可能有风险维度一大就要留神。最稳妥的做法是每轮相加后立即取模不要等全部累加完再统一取模。如果模数接近1e18甚至要用__int128过渡再转回long long。还有一个常用优化当矩阵维度固定且很小比如2×2或3×3可以把矩阵乘法展开成具体表达式省掉循环开销。这种写法很丑但在反复调用时会明显提速。3.3 用Python在本地验证建模是否正确C模板适合比赛但调试时我更推荐先用Python验证思路因为Python写起来快而且不会因为指针和数组越界干扰判断。下面这个简易版本足够用了def mat_mul(A, B, mod): n len(A) m len(B) p len(B[0]) C [[0] * p for _ in range(n)] for i in range(n): for k in range(m): if A[i][k] 0: continue for j in range(p): C[i][j] (C[i][j] A[i][k] * B[k][j]) % mod return C def mat_pow(A, e, mod): n len(A) R [[1 if i j else 0 for j in range(n)] for i in range(n)] while e: if e 1: R mat_mul(R, A, mod) A mat_mul(A, A, mod) e 1 return R我在带训练时给同学定了一个死规矩任何新题先用Python验证转移矩阵在n1到n20的结果再写C交比赛。为什么是20因为20以内的朴素递推可以手算或暴力循环验证一旦前20项全对矩阵建模基本就稳了。Python版本很慢但验证期要的不是速度是正确率等验证完再用C模板替换实际多花的时间不超过五分钟。4. 实战中踩过的坑与排查思路4.1 常见矩阵快速幂错误速查表矩阵快速幂的报错模式其实很集中训练课上我列过一张内部速查表这次整理出来症状可能原因排查与修复求出的F(n)变成F(n1)或F(n-1)快速幂的指数多乘或少乘一次打印n2、3时的结果确定初始向量和幂次对应关系矩阵乘法结果完全不对循环行列顺序写反矩阵写成了转置用2×2矩阵手算验证一遍再写代码取模后出现负数递推式里有减法取模没有调整运算后先加MOD再取模程序超时每次询问都重新pow或矩阵维度太大预处理幂、合并多询问维度大时尝试向量乘矩阵快速幂结果始终是零矩阵初始结果矩阵设置成全零结果矩阵应初始化为单位矩阵边界n0或n1出错没有单独处理特判先特判n 2再进快速幂带n²的递推结果错n、n²、常数缺少状态位把n、n²、1全部纳入状态向量重新推导矩阵模数很大时答案溢出乘法累加超出long long精度使用__int128或及时分段取模这些坑每一个我都见过其中最隐蔽的是状态向量方向的左右混淆。代码里矩阵数组a[i][j]表达第i行第j列但公式里“左乘”还是“右乘”并不是写代码时能看出来的必须在注释里固定一种写法列向量统一写成v_new M × v_old。如果一会儿写成行向量v_new v_old × M整个转移矩阵就要转置很容易错乱。4.2 典型案例复盘图上恰好走k步的路径数这次训练课上有一道经典题很能说明矩阵快速幂的建模过程。题目大意是给定一个n个点m条边的有向图求从点s到点t恰好走k条边的方案数k最大到1e18结果取模。n通常不超过50m也不算大。做法是直接用邻接矩阵MM[i][j]表示i到j的边数有边记1无记0。答案就是(M^k)[s][t]。为什么这么写看k2的情况(M²)[s][t]等于对每个中间点mid求和M[s][mid]×M[mid][t]这正是从s经过任意一个中间点再到t的两步路径方案数。k3时就再乘一次M相当于继续枚举第三个中间点每一步都由矩阵乘法的计数语义自动完成。用数学归纳法可以证明k次乘法的结果就表示恰好k步的路径数。需要注意两个细节。第一k0时矩阵的0次幂是单位矩阵所以从s到t恰好0步的方案数只有当st时为1这正好对应“原地不动”的唯一路径。第二如果图是无向的邻接矩阵本身是对称的但转移矩阵照样能用。这类题扩展也很多比如路径长度不超过k时可以加一个虚拟节点或把“停住”状态也建模进矩阵不过那是后面进阶的内容了。4.3 典型案例复盘求前n项和的递推另一道高频题是求斐波那契的前n项和S(n) F(1)F(2)...F(n)。最直接的方案是套恒等式S(n)F(n2)-1只用一次矩阵快速幂算F(n2)就完事。但训练课上我更推荐用增广状态的方法演示一遍因为它能帮你理解“把目标量当作状态之一”的通用思路。构造状态向量v(n) [F(n), F(n-1), S(n)]^T从v(n-1)推到v(n)F(n) F(n-1) F(n-2)F(n-1) F(n-1)S(n) S(n-1) F(n) S(n-1) F(n-1) F(n-2)所以转移矩阵M是[1, 1, 0] [1, 0, 0] [1, 1, 1]初始向量v(1) [1, 0, 1]^T因为F(1)1F(0)0S(1)1。验证一下M² × v(1)可以得到[2,1,4]^T即F(3)2、F(2)1、S(3)4手算完全一致。其实这道题用增广状态做没有恒等式快但它展示了一种很值钱的思维你需要的量没出现在递推式里那就把它作为独立维度放进去。前缀和、区间和、累计次数这类问题很多都能用同一招解决。4.4 一个降低调试成本的独家小习惯最后说一个我自己带训多年攒下的习惯任何转移矩阵第一次写完不要直接交评测先做一个“小n验证”。具体做法是用朴素写法跑出前10到20项再用矩阵快速幂跑同样范围逐项对拍。这一步能过滤掉绝大部分建模错误包括状态向量方向、幂次差一、矩阵某行列漏写等常见问题。我见过太多同学在矩阵快速幂上卡一整晚最后发现只是某一行系数写反。原因是他们把调试重心放在快速幂函数上而不是矩阵本身。实际上快速幂的代码是成熟模板不容易出问题真正需要反复推理的是“这个递推式到底该怎么变成矩阵”。统一约定也很重要。我们内部的标准写法是状态一律列向量转移一律“新状态 矩阵 × 旧状态”矩阵数组a[i][j]严格表示状态i如何由状态j贡献系数。团队协作时这个约定能避免“我看你的代码像转置你看我的代码像左乘”的混乱。我个人带完这期训练班最大的感受是矩阵快速幂是那种“建模对了代码半小时就完事”的题。真正花时间的不是函数怎么写而是递推式怎么塞进矩阵里。所以后来在训练里定了一条规矩专题训练的前两题必须手推转移矩阵写完先对拍对拍通过再去想优化。这个小习惯帮我们减少了大量调整样例的时间。如果后续有时间我可以把AC自动机上的矩阵DP和带环图上的计数这两类进阶题也写成单独复盘它们其实是这个专题的自然延伸。
阅读完成 · 觉得有帮助?
咨询建站