资讯详情

资讯详情

矩阵连乘问题:动态规划从入门到实战的完整推导

矩阵连乘这个问题我在刚接触算法设计与分析那会儿第一反应就是这不就是把一堆矩阵按顺序乘起来吗有什么好讲的结果真动手写了个朴素递归去跑才发现自己太天真了。当矩阵链长度稍微上去一点比如十几个矩阵程序直接卡死递归树膨胀得比想象中快得多。后来才明白矩阵连乘问题之所以被反复拿出来讲是因为它把动态规划最核心的几个思想——最优子结构、重叠子问题、状态转移、区间划分——全都浓缩在了一个非常直观的场景里。它不只是一个算法题更像是理解动态规划的一把钥匙。这篇文章我想从一个实际写代码的人的角度把矩阵连乘问题从头到尾拆一遍。包括它到底在解决什么问题、为什么朴素方法不行、动态规划的状态怎么定义、转移方程怎么推、代码怎么写、边界怎么处理、以及实际跑起来容易踩哪些坑。适合正在学算法设计与分析的学生也适合工作中偶尔需要处理矩阵运算优化、想重新捡起动态规划思路的开发者。我会尽量把每一步的“为什么”讲清楚而不是只丢一个公式出来。1. 矩阵连乘问题到底在解决什么1.1 从矩阵乘法的基本代价说起先把最基础的东西捋清楚。两个矩阵相乘比如 A 是 p×q 的矩阵B 是 q×r 的矩阵它们相乘得到 C 是 p×r 的矩阵。这个乘法能做的前提是 A 的列数等于 B 的行数也就是那个 q 要对上。计算 C 中每一个元素需要做 q 次乘法和 q 次加法而 C 一共有 p×r 个元素所以总的标量乘法次数是 p×q×r。这个 p×q×r 就是衡量矩阵乘法代价的核心指标。为什么只数乘法次数因为在计算机里乘法通常比加法慢而且在算法分析中乘法次数和加法次数是同阶的数一个就够了。所以后面我们讨论“代价”默认说的就是标量乘法次数。现在问题来了。假设有三个矩阵 A1、A2、A3维度分别是 10×100、100×5、5×50。如果按 ((A1A2)A3) 的顺序算先算 A1A2代价是 10×100×5 5000得到一个 10×5 的矩阵再和 A3 相乘代价是 10×5×50 2500。总共 7500 次乘法。如果换个顺序按 (A1(A2A3)) 算先算 A2A3代价是 100×5×50 25000得到一个 100×50 的矩阵再和 A1 相乘代价是 10×100×50 50000。总共 75000 次乘法。同样的三个矩阵同样的最终结果仅仅因为加括号的顺序不同计算量差了十倍。这就是矩阵连乘问题的核心给定一个矩阵链怎么加括号也就是怎么确定相乘顺序使得总的标量乘法次数最少。1.2 为什么加括号的顺序会影响结果这里有个关键点容易被忽略矩阵乘法虽然满足结合律也就是 (AB)C A(BC)最终结果是一样的但它不满足交换律你不能随便调换矩阵的先后位置。所以你能做的只是决定“先算哪一对相邻的”也就是在哪些位置加括号而不能改变矩阵链本身的顺序。那为什么加括号会影响代价本质原因在于两个矩阵相乘会产生一个中间结果矩阵这个中间矩阵的维度决定了后续乘法的代价。你先把哪两个乘起来就会先“消掉”一个维度剩下的中间矩阵形状不同后续每一步的代价就不同。用生活化的类比来说这就像你要搬一堆不同大小的箱子箱子的搬运顺序会影响你中途需要多大的临时堆放空间和多少趟搬运。先合并小的可能中间产物小、后续省力先合并大的中间产物巨大后面每一步都拖累。矩阵连乘要做的就是找到那个“最省力的合并顺序”。1.3 问题的形式化描述把上面的直觉形式化一下。给定 n 个矩阵 A1, A2, ..., An其中 Ai 的维度是 p(i-1) × p(i)。注意这里维度是用一个数组 p 来表示的p 有 n1 个元素因为相邻矩阵共享一个维度。比如 A1 是 p0×p1A2 是 p1×p2以此类推。这个用一维数组表示维度链的技巧非常关键它把 n 个矩阵的维度信息压缩成了 n1 个数后面写代码会一直用到。问题是确定一个完全加括号的方式也就是确定一个乘法顺序使得计算 A1A2...An 所需的总标量乘法次数最少。所谓完全加括号是指最终可以化归为两个子链的乘积而每个子链又递归地完全加括号。这里要强调一个概念最优子结构。如果整个链的最优加括号方案在某个位置 k 把链分成 A1...Ak 和 A(k1)...An 两部分那么这两部分各自的加括号方案也必须是它们各自的最优方案。为什么因为如果某一边存在更省的计算方式换上去只会让整体更省那原来的方案就不是最优的了。这个性质是动态规划能用的前提也是整个解法成立的基石。2. 朴素方法为什么不可行2.1 穷举所有加括号方案的数量最直接的想法是把所有可能的加括号方式都试一遍取代价最小的那个。那到底有多少种加括号方式这个数量其实是一个卡特兰数。n 个矩阵的完全加括号方案数是第 n-1 个卡特兰数记作 C(n-1)。卡特兰数的增长非常快。C(1)1C(2)2C(3)5C(4)14C(5)42C(10)16796C(15) 已经接近一千万C(20) 是几十亿的量级。也就是说当矩阵个数到 20 个左右穷举法就已经完全跑不动了。而实际工程中矩阵链长度几十上百都很常见穷举根本不可能。我当初写那个朴素递归就是隐式地在枚举所有加括号方案。代码大概是这样对区间 [i, j]枚举所有分割点 k递归计算左边和右边再加上合并代价。这个递归没有记忆化所以同一个子区间被反复计算了无数次复杂度是指数级的。跑十几个矩阵就卡死一点都不奇怪。2.2 重叠子问题的发现真正让动态规划登场的原因是这些子问题被大量重复计算。还是拿 [i, j] 这个区间来说无论它是作为哪个更大区间的左半部分还是右半部分被调用它本身的“最优加括号代价”是固定的只跟 i 和 j 有关跟外层怎么切无关。举个具体的例子。假设有 5 个矩阵我们在计算整个链 [1,5] 的时候会需要 [1,2]、[2,5]、[1,3]、[3,5] 等等子区间。而在计算 [1,4] 的时候又会需要 [1,2]、[2,4]、[1,3]、[3,4]。你会发现 [1,2]、[1,3] 这些子区间被反复用到。朴素递归每次都重新算一遍白白浪费了大量时间。这就是典型的重叠子问题。解决办法也很自然把算过的子区间结果存起来下次直接查表。这就是记忆化搜索也是动态规划自顶向下写法的核心。而如果进一步按区间长度从小到大填表就变成了自底向上的动态规划。两种写法本质一样只是方向不同。2.3 指数级到多项式级的跨越加上记忆化之后复杂度会发生质的变化。子问题的数量是多少所有可能的区间 [i, j]其中 1 ≤ i ≤ j ≤ n数量是 n(n1)/2也就是 O(n²)。每个子问题在求解时需要枚举分割点 k从 i 到 j-1最多 n 种选择。所以总的时间复杂度是 O(n³)。空间复杂度是 O(n²)因为要存所有子区间的结果。从指数级降到多项式级这就是动态规划在这个问题上的威力。n100 的时候n³ 是一百万随便跑而指数级的话宇宙毁灭都算不完。这个对比非常震撼也是为什么矩阵连乘被当作动态规划的经典入门案例。3. 动态规划解法的完整推导3.1 状态定义m[i][j] 的含义动态规划的第一步永远是定义状态。在矩阵连乘里我们定义 m[i][j] 表示计算子链 AiAi1...Aj 所需的最少标量乘法次数。注意这里的 i 和 j 是矩阵的下标从 1 开始。为什么这样定义因为整个问题的答案就是 m[1][n]。而且这个定义天然满足最优子结构m[i][j] 的最优值一定是由某个分割点 k 把区间分成 [i,k] 和 [k1,j] 两部分再加上这两部分合并的代价得到的。子问题 m[i][k] 和 m[k1][j] 都是同样形式的问题只是规模更小。这里有个细节要注意m[i][i] 表示只有一个矩阵的情况不需要做任何乘法所以 m[i][i] 0。这是递推的边界条件也是填表的起点。3.2 状态转移方程的推导过程现在推导转移方程。对于区间 [i, j]假设最优分割点在 ki ≤ k j那么计算过程是先算左半部分 Ai...Ak代价是 m[i][k]再算右半部分 A(k1)...Aj代价是 m[k1][j]最后把这两个结果矩阵相乘。左半部分的结果矩阵维度是多少Ai...Ak 乘完之后得到的是一个 p(i-1) × p(k) 的矩阵。右半部分 A(k1)...Aj 乘完得到的是 p(k) × p(j) 的矩阵。这两个矩阵相乘代价是 p(i-1) × p(k) × p(j)。所以对于固定的 k总代价是m[i][k] m[k1][j] p(i-1) × p(k) × p(j)而我们要在所有可能的 k 中取最小值于是转移方程是m[i][j] min{ m[i][k] m[k1][j] p(i-1)×p(k)×p(j) }其中 i ≤ k j这个方程是整个解法的核心。我第一次看的时候卡在那个 p(i-1)×p(k)×p(j) 上没想明白为什么是这三个数相乘。后来画了个图才清楚左半结果的行数是 p(i-1)列数是 p(k)右半结果的行数是 p(k)列数是 p(j)两者相乘中间维度 p(k) 对上代价就是 p(i-1)×p(k)×p(j)。想通这一点整个方程就顺了。3.3 填表顺序与区间长度递增有了转移方程接下来要确定计算顺序。m[i][j] 依赖于 m[i][k] 和 m[k1][j]这两个子区间的长度都比 [i,j] 短。所以我们必须先算短的区间再算长的区间。这就是按区间长度递增填表。具体来说先初始化所有 m[i][i] 0这是长度为 1 的区间。然后计算长度为 2 的区间 m[i][i1]它们只依赖长度为 1 的子区间。接着长度 3、4一直到长度 n。最后 m[1][n] 就是答案。这个填表顺序非常关键。如果顺序搞错比如按 i 从小到大、j 从小到大填可能会用到还没算出来的子区间结果就错了。我见过不少人在这里翻车代码逻辑看着没问题但结果就是不对最后发现是循环顺序写反了。除了记录最小代价通常还需要一个辅助表 s[i][j] 来记录最优分割点 k。这样不仅能算出最小代价还能回溯出具体的加括号方案。s[i][j] 存的就是让 m[i][j] 取到最小值那个 k。4. 代码实现与关键细节4.1 自底向上动态规划的完整实现下面给出一个完整的 Python 实现。我用的是 1-based 下标来对应矩阵编号但 Python 列表是 0-based所以代码里会做一些偏移处理这点要特别注意。def matrix_chain_order(p): p: 维度数组长度为 n1表示 n 个矩阵的维度链 返回: (m, s)m 是最小代价表s 是最优分割点表 n len(p) - 1 # 矩阵个数 # m[i][j] 表示计算 Ai...Aj 的最少乘法次数 # 使用 (n1) x (n1) 的表下标 1..n 有效 m [[0] * (n 1) for _ in range(n 1)] s [[0] * (n 1) for _ in range(n 1)] # l 是区间长度从 2 开始 for l in range(2, n 1): for i in range(1, n - l 2): j i l - 1 m[i][j] float(inf) for k in range(i, j): cost m[i][k] m[k 1][j] p[i - 1] * p[k] * p[j] if cost m[i][j]: m[i][j] cost s[i][j] k return m, s这段代码有几个地方值得展开说。外层循环 l 是区间长度从 2 到 n。为什么从 2 开始因为长度 1 的区间 m[i][i] 已经在初始化时设为 0 了不需要计算。内层 i 的范围是 1 到 n-l1这样 j il-1 不会越界。最内层 k 从 i 到 j-1枚举所有分割点。那个 cost 的计算就是转移方程的直接翻译。p[i-1] * p[k] * p[j] 就是合并代价。注意这里的下标p[i-1] 是左半结果的行数p[k] 是左半结果的列数也是右半结果的行数p[j] 是右半结果的列数。这个对应关系一定要搞清楚写错了结果就全乱了。4.2 回溯构造最优加括号方案光算出最小代价还不够通常我们还想知道具体怎么加括号。这时候 s 表就派上用场了。s[i][j] 记录了区间 [i,j] 的最优分割点我们可以递归地构造出加括号的字符串。def print_optimal_parens(s, i, j): if i j: return fA{i} else: k s[i][j] left print_optimal_parens(s, i, k) right print_optimal_parens(s, k 1, j) return f({left} x {right})这个递归很直观如果区间只有一个矩阵直接返回矩阵名否则取出分割点 k递归处理左右两边然后用括号包起来。最终得到的字符串就是最优的加括号方案。我实测过对于 6 个矩阵的经典例子维度链是 [30, 35, 15, 5, 10, 20, 25]算出来的最小代价是 15125最优方案是 ((A1 x (A2 x A3)) x ((A4 x A5) x A6))。你可以自己跑一遍验证这个例子在很多教材里都出现过是个很好的对照。4.3 记忆化搜索版本对比除了自底向上还可以写成自顶向下的记忆化搜索。两种写法思路一样但代码风格不同适合不同习惯的人。def matrix_chain_memo(p): n len(p) - 1 m [[-1] * (n 1) for _ in range(n 1)] s [[0] * (n 1) for _ in range(n 1)] def lookup(i, j): if m[i][j] 0: return m[i][j] if i j: m[i][j] 0 return 0 m[i][j] float(inf) for k in range(i, j): cost lookup(i, k) lookup(k 1, j) p[i - 1] * p[k] * p[j] if cost m[i][j]: m[i][j] cost s[i][j] k return m[i][j] lookup(1, n) return m, s记忆化版本的好处是代码更贴近转移方程的数学形式读起来更自然而且只计算真正用到的子问题。缺点是递归有栈深度限制n 特别大的时候可能爆栈。自底向上版本没有这个问题但会计算所有子问题哪怕有些用不到。实际选择看场景一般 n 在几百以内两种都行。对比维度自底向上 DP记忆化搜索代码风格循环填表递归加缓存计算范围所有子问题仅用到的子问题栈风险无n 大时可能爆栈可读性循环顺序需理解贴近数学公式适用场景n 较大、追求稳定n 较小、追求直观5. 实操中容易踩的坑与排查5.1 下标偏移导致的越界和错值这是最常见的坑没有之一。矩阵连乘的数学描述用的是 1-based 下标矩阵 Ai 的维度是 p(i-1) × p(i)。但 Python、Java、C 这些语言的数组是 0-based。如果你直接把数学公式翻译成代码而不做偏移处理很容易出现 p[-1] 这种越界或者取到错误的维度值。我的建议是要么像上面代码那样表开成 (n1)×(n1)下标 1 到 n 有效0 号位置空着不用要么全程用 0-based但转移方程里的下标要相应调整。两种都行关键是统一别混着用。我见过有人在 m 表上用 1-based在 p 数组上用 0-based结果 p[i-1] 和 p[i] 来回搞混调了半天才发现。5.2 循环顺序写反的典型症状前面提过填表要按区间长度递增。如果写成按 i 递增、j 递增的双重循环会怎样比如for i in range(1, n 1): for j in range(i 1, n 1): # 计算 m[i][j]这样写的问题在于计算 m[i][j] 时需要的 m[k1][j] 可能还没算出来。因为 k1 可能大于 i而外层 i 是从小到大当 i 较小时较大的 i 对应的行还没处理到。结果就是用到初始值 0算出来的代价偏小最终答案错误。症状是程序不报错能跑出结果但结果比正确答案小。这种“静默错误”最难查。排查方法很简单拿一个小例子手算一遍对比程序输出。如果对不上先检查循环顺序。5.3 代价溢出与数据类型选择p(i-1) × p(k) × p(j) 这个乘积当矩阵维度很大时可能超出 32 位整数的范围。比如维度都是几千三个数相乘就是几十亿超过 int 上限。在 C/C 里要用 long longJava 里要用 longPython 倒是不用担心整数自动扩展。另外m[i][j] 的初始值要设成一个足够大的数代表无穷大。用 float(inf) 最稳妥或者用一个肯定大于任何可能代价的数比如 sum(p) 的三次方之类的。但别用 999999 这种拍脑袋的数万一真实代价超过它结果就错了。5.4 常见问题速查表问题现象可能原因排查方向结果比正确答案小循环顺序错误子问题未算先用检查是否按区间长度递增数组越界报错下标偏移处理不当统一 1-based 或 0-based结果完全不对合并代价公式写错核对 p[i-1]*p[k]*p[j]大 n 时结果异常整数溢出换 long long 或大整数回溯方案括号错位s 表记录时机不对确认在更新最小值时同步更新 s程序卡死用了无记忆化的朴素递归加缓存或改自底向上提示调试矩阵连乘问题时先用 3 到 4 个矩阵的小例子手算验证确认无误后再上大规模数据。小例子能暴露绝大多数逻辑错误。6. 从矩阵连乘看动态规划的通用套路6.1 区间 DP 的识别特征矩阵连乘是区间动态规划的典型代表。所谓区间 DP就是状态定义在一个区间 [i,j] 上转移时枚举区间内的分割点把大区间拆成两个小区间。识别这类问题的特征很明显问题可以分解成两个子问题的组合子问题对应原问题的连续子段且最优解由子问题最优解组合而成。一旦识别出是区间 DP套路基本固定定义 m[i][j] 为区间最优值枚举分割点 k转移方程形如 m[i][j] opt{ m[i][k] m[k1][j] 合并代价 }按区间长度递增填表。矩阵连乘、石子合并、多边形三角剖分、最长回文子序列都是这个套路。掌握一个其他的触类旁通。6.2 状态设计的通用思路动态规划最难的一步永远是状态设计。矩阵连乘的状态 m[i][j] 之所以好想是因为问题本身就带着“区间”这个结构。但在很多问题里区间结构不是显式的需要你自己去发现。一个通用的思路是问自己“子问题是什么”。在矩阵连乘里子问题是“某一段连续矩阵链的最优代价”。这个“连续”很关键它保证了子问题可以用两个下标 [i,j] 来描述。如果子问题不是连续的那可能就需要更复杂的状态比如集合状压或者树形结构。另一个思路是看“决策是什么”。矩阵连乘的决策是“在哪里切一刀”。决策的枚举范围决定了转移的复杂度。这里决策是 O(n) 种所以总复杂度 O(n³)。如果决策是 O(1) 种那复杂度就降到 O(n²)。分析决策空间是估计算法复杂度的关键。6.3 复杂度优化的可能方向O(n³) 的矩阵连乘在 n 到几百的时候还能接受但 n 上千就有点吃力了。有没有更快的算法其实是有的。有一个叫 Hu-Shing 的算法能把矩阵连乘优化到 O(n log n)但实现相当复杂而且常数很大实际中很少用。对于绝大多数场景O(n³) 已经够用。如果只是想在实际项目中减少矩阵乘法的代价与其纠结最优加括号不如考虑别的优化比如用更快的矩阵乘法库、利用稀疏性、或者干脆改变算法整体结构。矩阵连乘的最优加括号更多是理论意义和教学意义实际工程中矩阵链长度往往没那么夸张或者瓶颈根本不在加括号顺序上。注意不要为了追求理论最优而过度工程化。先测量实际瓶颈再决定要不要上复杂算法。我见过有人花大力气实现了一个复杂优化结果发现瓶颈在数据读取上白忙一场。6.4 一个容易忽略的实操心得最后分享一个我在实际写这类 DP 时的小习惯把转移方程先用注释写在代码里再逐行翻译成代码。比如# m[i][j] min over k in [i, j-1] of # m[i][k] m[k1][j] p[i-1]*p[k]*p[j]这样写的好处是代码和数学公式一一对应检查的时候直接对照不容易漏项或写错下标。尤其是合并代价那一项很容易漏掉或者写错乘数。先写注释再写代码能省下大量调试时间。另外我习惯在填表完成后把 m 表和 s 表打印出来看看。对于小例子肉眼扫一遍就能发现异常比如某些位置明显偏大或偏小。这个习惯帮我抓出过好几次循环边界的问题。表格可视化是调试 DP 的利器强烈建议养成。矩阵连乘这个问题表面上看是在算一个最优加括号顺序实际上是在训练一种思维方式如何把一个大问题拆成可复用的子问题如何定义状态让子问题之间有清晰的依赖关系如何按正确的顺序求解。这套思维方式才是动态规划真正值钱的地方。把矩阵连乘吃透再去看其他 DP 问题会发现很多都是换了个壳的同类问题。
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →