矩阵连乘问题Matrix Chain MultiplicationMCM是我在准备算法面试时反复遇到的一个经典动态规划题目。初次看它时我以为这不就是把矩阵乘起来嘛有什么可优化的真正缕清之后才发现这个看似简单的问题背后藏着一整套动态规划的标准思路、边界处理技巧和代码实现细节。今天这篇笔记我就把完整的推导过程、代码模板、常见坑点和应用场景一次讲清楚适合正在学数据结构与算法、准备机试或想彻底搞懂区间DP的读者。先说清楚矩阵连乘问题到底在解决什么给定n个矩阵的链我们要找到一种加括号的方式使得计算整个矩阵乘积时标量乘法的总次数最少。注意这里不要求真正算出乘积结果只求最少乘法次数和对应的括号方案。这个问题最吸引我的地方在于它看起来像是一个“顺序选择”问题实际却是一个“区间划分”问题能用动态规划优雅解决而且在深度学习、计算机视觉、图形学变换矩阵组合中有直接应用。1. 问题建模与暴力解法的真相1.1 矩阵乘法的代价是怎么算的要理解这个问题先得回到矩阵乘法本身。假设矩阵A是m×n的矩阵B是n×p的那A×B的结果矩阵是m×p的。为了算出结果里每一个元素我们要做n次乘法而结果矩阵里有m×p个元素所以总共需要的标量乘法次数是m×n×p。举个例子A是2×3的B是3×4的那么A×B的乘法次数就是2×3×424次。这个代价公式是整个矩阵连乘问题的基石后面的所有推导都建立在这条公式上。三条矩阵相乘时情况开始变得有趣了。假设有三个矩阵A1是10×100的A2是100×5的A3是5×50的。按(A1×A2)×A3的顺序计算先算A1×A2需要10×100×55000次乘法得到的结果矩阵是10×5的再乘以A3需要10×5×502500次乘法总共7500次。如果按A1×(A2×A3)的顺序先算A2×A3需要100×5×5025000次乘法得到10×100的中间结果再乘以A1需要10×100×5050000次总共75000次。同样三个矩阵只是因为加括号的顺序不同乘法次数相差整整10倍。这个例子让我第一次意识到矩阵乘法虽然满足结合律但计算代价对运算顺序极其敏感。矩阵连乘问题要做的就是在所有可能的加括号方式里找到代价最低的那一种。1.2 为什么暴力搜索不可行既然要找最优括号方案最直接的想法就是枚举所有可能的加括号方式分别计算代价再取最小值。但这里藏着一个指数爆炸的问题。n个矩阵相乘所有可能的加括号方式数量是一个卡特兰数Catalan number记作C(n-1)。它的递推关系是C(1)1C(n)ΣC(k)×C(n-k)k从1到n-1。前几项是C(1)1C(2)1C(3)2C(4)5C(5)14C(6)42C(7)132C(8)429C(9)1430。看着增长不算快但当n20时C(19)已经超过10亿。也就是说当矩阵数量到二三十个时暴力枚举已经完全没有可能在合理时间内完成。这里要回应一个我经常看到的误解很多人以为暴力解法就是全排列矩阵先做全排列再乘。实际上括号方案和排列没有直接关系矩阵的顺序是固定的我们只能决定在哪里加括号这个数量就是卡特兰数。全排列n!会比卡特兰数更夸张但即便用卡特兰数来算也已经足够让暴力解法在n稍大时彻底失效。1.3 贪心思路为什么第一个被排除在接触动态规划之前我第一反应是贪心是不是每次都选择相邻矩阵中乘法代价最小的一对先乘这样每一步都局部最优最后全局是不是也最优实测下来完全不对。还是用刚才那组数据A110×100A2100×5A35×50。按贪心思路先看A1×A2的代价是5000A2×A3的代价是25000选择先乘A1×A2得到中间矩阵是10×5的再乘A3代价是10×5×502500总代价7500。这个结果碰巧是好的。但如果换一组数据A1100×1A21×100A3100×1。贪心先找A1×A2代价是100×1×10010000A2×A3也是100×1×10010000选哪个都一样。假设选A1×A2得到中间矩阵100×100再乘A3代价是100×100×110000总代价20000。但如果我们先乘A2×A3代价10000得到中间1×1的矩阵再乘A1的代价是100×1×1100总代价10100。贪心直接多算了一倍的乘法次数。这个反例说明了一个核心问题矩阵连乘本身的乘数之间不是独立的先乘哪一对会影响后续中间矩阵的维度而中间矩阵的维度又直接决定下一步的代价。全局最优解要求我们综合考虑所有划分可能这正是动态规划能发挥作用的地方。2. 动态规划核心思路状态与转移方程2.1 状态定义dp[i][j]的含义既然要避免重复计算动态规划的第一步是定义状态。这个问题的经典状态定义是用dp[i][j]表示从第i个矩阵连乘到第j个矩阵即Ai×Ai1×...×Aj所需的最少标量乘法次数。这个定义是典型的“区间DP”思路。我们不去想整个矩阵链怎么乘而是把它拆成任意一个连续的区间[i, j]先求出这个区间的最优代价再组合成更大的区间。为什么状态必须是连续的区间因为矩阵乘法的结合性决定了任何一次括号划分最终都把原问题分成两部分左边是连续的若干矩阵右边也是连续的若干矩阵。整个链的划分树展开后每一棵子树都对应一个连续的矩阵区间。所以只要把所有连续区间的最优解都算出来最后dp[1][n]就是整个问题的最优解。这里还要明确下标边界。假设有n个矩阵A1, A2, ..., An我们用数组p来存储维度信息p[i-1]表示Ai的行数p[i]表示Ai的列数。也就是说矩阵Ai的维度是p[i-1]×p[i]。这样一来相邻矩阵的维度天然衔接Ai的列数p[i]等于Ai1的行数p[i-1]的后一个下标。p的长度是n1下标从0到n。我个人觉得下标对应关系是初学矩阵连乘最容易绕晕的地方很多代码bug都出在这里。可以用一句话记住p数组是n1个元素的维度序列第i个矩阵的行和列由p[i-1]和p[i]决定。2.2 状态转移方程怎么推有了dp[i][j]的定义下一步就是推导状态转移方程。考虑区间[i, j]我们在这个区间内选择一个分割点k把矩阵链分成两部分左边的[i, k]和右边的[k1, j]分别算出最优代价再加上把左右两部分的结果矩阵相乘的代价就得到以k为分割点时这个区间的总代价。合并代价怎么算左边部分[i, k]乘法完成后得到一个维度为p[i-1]×p[k]的矩阵右边部分[k1, j]乘法完成后得到一个维度为p[k]×p[j]的矩阵两者相乘的标量乘法次数是p[i-1]×p[k]×p[j]。所以整个区间的总代价是dp[i][k] dp[k1][j] p[i-1]×p[k]×p[j]我们要在k从i到j-1的所有可能取值里找最小值所以状态转移方程是dp[i][j] min{dp[i][k] dp[k1][j] p[i-1]×p[k]×p[j]}i ≤ k j这个方程写出来很简单但背后的含义值得反复咀嚼。它把一个大问题拆分成两个规模更小的子问题子问题互相独立没有重叠依赖的循环因此满足动态规划的最优子结构性质。2.3 为什么循环顺序要对区间长度从小到大状态转移方程确定后最难理解的部分来了循环顺序。很多人会自然地想既然dp[i][j]依赖dp[i][k]和dp[k1][j]那我用三重循环i、j、k遍历不就行了吗实测下来不行因为dp[i][j]依赖的区间长度比[i, j]本身短如果外层循环从小到大枚举i内层循环枚举j可能遇到dp[i][k]或者dp[k1][j]还没计算的情况。正确的做法是外层循环枚举区间长度len从2开始到n中层循环枚举区间起点i那么区间终点j ilen-1内层循环枚举分割点k从i到j-1。这样保证计算任何一个长度为len的区间时所有长度小于len的区间都已经提前算好了依赖关系永远是从短区间指向长区间。这里的直观理解可以类比为填一张二维表格我们是从左上角往右下角对角线方向一格一格地填而不是一行一行地填。因为dp[i][j]依赖的是表格中位于它左侧和下方的两组格子按区间长度递增的顺序填充才能保证这些格子已经就绪。我建议初学者动手把表格画出来把每个格子依赖的格子标出来循环顺序就一目了然了。这一步想通了区间DP基本就入门了。3. 代码实现与输出最优括号方案3.1 核心代码计算最小乘法次数先看最核心的计算部分。我习惯用C写算法题这里给出一个可以直接跑的完整函数片段注释写得比较细#include bits/stdc.h using namespace std; const int INF 0x3f3f3f3f; // p[0..n]p[i-1]是第i个矩阵的行数p[i]是列数 // n矩阵个数 int matrixChainOrder(vectorint p, int n, vectorvectorint s) { // dp[i][j]第i到第j个矩阵的最小乘法次数 vectorvectorint dp(n 1, vectorint(n 1, 0)); s.assign(n 1, vectorint(n 1, 0)); // s[i][j]记录分割点 for (int len 2; len n; len) { // len区间长度 for (int i 1; i len - 1 n; i) { // i区间起点 int j i len - 1; // j区间终点 dp[i][j] INF; for (int k i; k j; k) { // k分割点 int cost dp[i][k] dp[k1][j] p[i-1] * p[k] * p[j]; if (cost dp[i][j]) { dp[i][j] cost; s[i][j] k; // 记录最优分割点便于回溯 } } } } return dp[1][n]; }关键点在于p[i-1] * p[k] * p[j]的写法。很多初学者会写成p[i] * p[k] * p[j1]之类这往往是因为维度下标没理清楚。只要记住左边区间的结果矩阵是p[i-1]×p[k]的右边区间的结果矩阵是p[k]×p[j]的合并代价就是三者相乘这个式子就不会错。3.2 回溯分割点输出加括号方案光算出最小乘法次数还不够实际应用中往往需要输出具体的加括号方案。这需要在状态转移时额外记录每个区间的最优分割点s[i][j]然后通过递归回溯输出。回溯的思路是对于区间[i, j]如果i等于j说明只有一个矩阵直接输出这个矩阵的名字否则先输出左括号递归输出[i, s[i][j]]和[s[i][j]1, j]再输出右括号。代码如下void printOptimal(int i, int j, vectorvectorint s, vectorstring names) { if (i j) { cout names[i]; return; } cout (; printOptimal(i, s[i][j], s, names); printOptimal(s[i][j] 1, j, s, names); cout ); }这里的names可以是一个字符串数组比如{, A1, A2, ..., An}让输出更直观。回溯时要注意递归边界必须是ij而不是ij因为区间长度至少为1。我实际测试过一个经典数据矩阵维度分别是30×35、35×15、15×5、5×10、10×20、20×25也就是p数组为{30,35,15,5,10,20,25}n6。程序输出的最优括号方案是((A1(A2A3))((A4A5)A6))最小乘法次数是15125次。这个结果和经典算法教材《算法导论》里的例子完全一致可以作为自测用例。3.3 完整测试与运行效果演示写代码时我习惯顺手做一个自测确保算法正确后再拿去应对更复杂的数据。这里给出一段简单的完整测试代码int main() { vectorint p {30, 35, 15, 5, 10, 20, 25}; int n p.size() - 1; vectorvectorint s; int ans matrixChainOrder(p, n, s); cout 最小乘法次数: ans endl; vectorstring names {, A1, A2, A3, A4, A5, A6}; printOptimal(1, n, s, names); cout endl; return 0; }运行结果最小乘法次数: 15125 ((A1(A2A3))((A4A5)A6))这组数据虽然规模小但很能说明问题。如果按照最直观的思路从左往右不加括号地硬算也就是((((A1A2)A3)A4)A5)A6乘法次数是多少呢我算下来是30×35×15 30×15×5 30×5×10 30×10×20 30×20×25 15750 2250 1500 6000 15000 40500次。对比最优解15125次差距超过两倍半。这就直观说明了矩阵连乘优化的威力。4. 空间与复杂度分析为什么能更快4.1 时间复杂度与空间复杂度的详细评估矩阵连乘动态规划解法的时间复杂度是O(n^3)空间复杂度是O(n^2)。这个复杂度具体怎么来的三重循环外层枚举区间长度内层枚举起点最内层枚举分割点三个维度都是n的量级相乘就是n的三次方。不过实际计算次数可以更精确一点。设总计算量为Σ(len从2到n)Σ(i从1到n-len1)(len-1)。这个求和算下来大约是n的三次方除以6。也就是说当n100时大约执行16.7万次状态更新当n1000时约1.67亿次。这个量级在竞赛或者工程场景里用C或者Java处理n500左右的矩阵链是很轻松的。空间方面dp表和s表各需要(n1)×(n1)的二维数组对于n1000来说大约是8MB如果每个元素用int实际上dp和s各约4MB共8MB在一般的算法环境里完全够用。但要注意如果使用long long类型空间会翻倍需要提前规划。4.2 空间优化的可能性与局限性有些教程提到可以把空间从O(n^2)优化到O(n)但矩阵连乘的区间DP不太容易做到这一点。原因在于dp[i][j]的更新依赖的是任意两个子区间的组合不像某些序列DP只依赖前一个状态很难只保留一维数组完成全部更新。不过有个折中方案如果只需要最小乘法次数而不需要回溯括号方案可以把s表省掉只保留dp表空间从两个二维数组减到一个能省约一半内存。这在n特别大的时候还是很有用的。另外可以用备忘录法记忆化搜索实现同样的时间复杂度但实际工程中可能存在递归栈溢出的风险。我的建议是n小于500用自底向上DPn大于500但需要快速实现时考虑备忘录法配合显式栈或者直接上自底向上并用long long防溢出。4.3 经典优化技巧四边形不等式如果遇到n特别大的场景还有一个值得一提的优化四边形不等式优化。这个优化适用于满足特定性质的区间DP可以把时间复杂度从O(n^3)降到O(n^2)。四边形不等式的要求是对于i ≤ i ≤ j ≤ j满足dp[i][j] dp[i][j] ≤ dp[i][j] dp[i][j]。在矩阵连乘问题中这个性质是成立的因为代价函数满足四边形不等式的条件。利用这个性质每次枚举分割点k时范围可以从[i, j-1]缩小到[s[i][j-1], s[i1][j]]从而大幅减少内层循环次数。这个优化在竞赛题里偶尔用得着工程上如果矩阵数量特别大比如几千个矩阵可以考虑。但要注意它的实现比普通版本复杂要额外维护s数组并且注意边界。我的建议是先把基础版写熟再学这个优化否则容易出错。5. 常见问题与错误排查那些我踩过的坑5.1 下标溢出与维度对应错误矩阵连乘最经典的bug就是维度下标写错。最常见的情况是把p[i-1] * p[k] * p[j]写成p[i] * p[k] * p[j1]或者把dp[k1][j]写成dp[k][j]。这类错误的特点是数据量小的时候可能碰巧输出正确一旦换一组数据就崩。排查方法很简单我建议写一个暴力枚举函数用递归枚举所有括号方案算出真实最优值然后和DP的结果做对拍。只要随机生成长度不超过8的矩阵维度数据对拍几次就能把下标问题暴露出来。这里提供一个对拍用的暴力函数// 暴力递归返回[i, j]区间的最优代价 long long bruteForce(vectorint p, int i, int j) { if (i j) return 0; long long best LLONG_MAX; for (int k i; k j; k) { best min(best, bruteForce(p, i, k) bruteForce(p, k 1, j) 1LL * p[i-1] * p[k] * p[j]); } return best; }把dp版本的结果和暴力版本的结果对比如果n6、7时已经出现不一致那一定是某个细节写错了。5.2 INF设置太小导致状态未更新很多人在初始化dp[i][j]时习惯用INT_MAX这本身没问题但如果后续做加法时出现整形溢出问题就来了。dp[i][k]加上dp[k1][j]再加上p[i-1]*p[k]*p[j]如果p的某个维度特别大比如1000×1000三个维度相乘就是10亿几个10亿再加起来很容易就超过INT_MAX导致加法溢出变成负数最终得到错误的最小值。我的解决方法是在涉及乘法或累加的地方统一用long long或者提前把INF设置成足够大的数比如0x3f3f3f3f约10.7亿然后所有dp数组都用long long存储。矩阵连乘的答案在最坏情况下规模有多大n个维度全是1000的矩阵乘法次数大约是1000×1000×1000×nn100时就是1000亿远远超过int能表示的范围。所以用long long是必须的这一点千万不能省。5.3 输入边界情况矩阵数量为1或2当n1时不需要任何乘法dp[1][1]0程序应该直接输出0并退出。当n2时只有一种括号方案就是直接乘dp[1][2]p[0]*p[1]*p[2]不需要进入内层k循环。这两种边界情况都要在代码里正确处理。有些同学在写循环时会把len直接从1开始枚举这样会导致dp[i][i]被重新赋值可能干扰后续状态计算。我的建议是len从2开始并且额外初始化dp[i][i]0通过vector默认值已经实现这样最保险。如果输入数据存在非法情况比如某个矩阵的行列数为0或负数应该在读入时直接判断并报错。矩阵连乘算法本身假设所有矩阵都是可乘的即相邻矩阵的维度必须匹配否则问题无解。这一点在实际工程中要注意。5.4 常见问题速查表现象可能原因解决方法小数据正常、大数据答案错误整数溢出所有dp和乘法改为long long输出结果比真实值小INF设置过大/过小或加法溢出为负INF设为0x3f3f3f3f用long long输出全是INFlen循环起点错误dp[i][j]未更新len从2开始确保k从i到j-1回溯输出括号异常s数组未正确记录分割点检查s[i][j]k是否在更新处执行大n时性能差未用四边形不等式优化或递归过深自底向上DP必要时四边形不等式在某些编译器上报错vector下标越界检查ilen-1 ≤ np下标不能超n这张表是我从实际调试经验中总结出来的遇到问题可以先对号入座。6. 实战应用矩阵连乘不只在算法题里6.1 图形学与计算机视觉中的多矩阵变换热词里提到“绕任意轴旋转后坐标形式七矩阵连乘”这正是矩阵连乘的一个典型应用场景。在计算机图形学和计算机视觉中一个物体的位姿变化往往由多个基本变换矩阵组合而成比如平移、绕X轴旋转、绕Y轴旋转、绕Z轴旋转、缩放等。如果我们要对大量顶点执行同一组变换操作这组变换矩阵的乘积顺序就决定了效率。举个具体的例子要对一个三维模型绕任意轴旋转经典的Rodrigues旋转公式可以分解为若干个基础矩阵的乘积。如果先组合好这组矩阵七矩阵连乘再对每个顶点执行一次矩阵乘法比每步变换都单独乘一遍要快得多。这时候矩阵连乘算法可以帮助我们找到最优的矩阵组合顺序减少组合阶段的计算量。在深度学习领域卷积神经网络中张量操作如转置、reshape、矩阵乘法也可以通过类似的矩阵链优化来减少计算开销。计算机视觉中的多视角几何、姿态估计等任务常常涉及多个矩阵变换的组合用矩阵连乘思路做优化能显著提升效率。6.2 动态规划与其他算法的对比选择矩阵连乘问题除了用动态规划还有别的解法路径比如备忘录法自顶向下递归记忆化、分治思想、贪心算法、甚至最终的Strassen矩阵乘法与其他快速矩阵乘法的结合。用表格对比一下这些方案方法时间复杂度空间复杂度适用场景优点缺点暴力枚举括号O(C(n-1)) ≈ O(4^n/n^1.5)O(1)n≤10实现简单指数爆炸递归分治指数级O(n)理解思路直观大量重复计算备忘录法O(n^3)O(n^2)只算需要的区间避免无效计算递归栈风险动态规划自底向上O(n^3)O(n^2)通用推荐稳定、可回溯方案空间占用较大四边形不等式优化O(n^2)O(n^2)大规模n省时间实现复杂贪心O(n^2)或O(n)O(1)不可行快结果错误仅作反例实战中我优先推荐自底向上的动态规划代码稳定、思路清晰、可回溯方案对绝大多数应用场景来说性能足够。如果以后遇到矩阵数量特别大的工业级场景再考虑四边形不等式优化。6.3 从一道题到一类题区间DP的通用模板矩阵连乘是区间DP的经典入门题掌握它之后很多类似问题都可以套用同一个框架比如石子合并、多边形三角剖分、括号匹配等。这类题目的共同特征是问题可以划分为连续区间的子问题状态转移时枚举分割点最终求全局最优。我的经验是遇到一个新问题先问三个问题状态能不能定义为区间区间合并时的代价函数是什么是否需要记录分割点用于回溯如果能回答清楚基本就可以套用区间DP模板了。这个模板的通用结构是初始化所有dp[i][i]0 for len 2 to n: for i 1 to n-len1: j i len - 1 dp[i][j] INF for k i to j-1: dp[i][j] min(dp[i][j], dp[i][k] dp[k1][j] cost(i,k,j))一旦掌握了这个模板学其他区间DP题会顺很多。矩阵连乘就像一把钥匙打开了区间动态规划的大门。我自己在刷题和写工程代码的过程中越来越觉得算法题的真正价值不在于背模板而在于理解问题本质后能灵活迁移。矩阵连乘的问题本质是“如何通过调整计算顺序来降低组合代价”这个思想在很多领域都有体现。当你某天在优化一组矩阵变换、调整一组流水线任务顺序甚至安排一次多阶段计算的执行计划时可能就会想起这个经典问题的解法。