量化交易中文教材

第 15a 章 动态规划(上):原理与经典问题

动态规划(dynamic programming, DP)是量化交易里最常用的算法思想之一:最优执行的拆单、带交易成本的调仓路径、美式期权在二叉树上的倒推定价、隐马尔可夫模型的状态解码,本质上都是动态规划。本章分上下两篇。上篇按原书顺序把方法讲透:四个设计步骤、两个适用要素、两种实现方式,以及钢条切割、矩阵链乘法、最长公共子序列、最优二叉搜索树四个经典问题。下篇(15b)专讲量化交易中的动态规划建模。

学习目标

  1. 说清动态规划与分治法的区别,熟练使用"刻画结构—递归定义—自底向上计算—构造解"四步法。
  2. 会用"剪切—粘贴"论证证明最优子结构,能判断子问题是否独立,并能举出最优子结构不成立的反例(最长简单路径)。
  3. 理解重叠子问题如何让朴素递归变成指数时间,掌握带备忘的自顶向下与自底向上两种实现,知道各自的适用场合。
  4. 能用"子问题数 × 每个子问题的选择数"或子问题图估算 DP 的运行时间。
  5. 独立写出钢条切割、矩阵链乘法、LCS、最优 BST 的 DP,并能重构最优解。
  6. 能用矩阵链乘法的思想为因子风险模型的矩阵计算选择最省的计算顺序。

读前导读

这一章在解决什么问题。 动态规划只有两句话:把大问题拆成小问题,而且这些小问题会反复出现;每个小问题第一次算出来就记下答案,以后再遇到直接查。 名字里的 "programming" 是"填表"的意思,跟写程序无关。

你其实在 CFA 里已经做过动态规划,只是没这么叫:用二叉树给美式期权定价。一棵 3 步的二叉树,到期日有 4 个结点,每个结点的价值直接由行权收益给出;然后往回推一步,每个结点的价值 = max(立即行权的价值, 两个子结点价值的贴现期望);再往回推……直到根。这里有两个关键点。第一,"先涨后跌"和"先跌后涨"到达的是同一个结点(树是重合的),这个结点的价值只算一次,两条路径共用——这就是"重叠子问题"。如果不共用、把每条路径分开算,\(n\) 步就要算 \(2^n\) 条路径;共用之后只有约 \(n^2/2\) 个结点。第二,每个结点的最优决策(行权还是继续持有)只取决于从这个结点往后的情况,与你是怎么走到这个结点的无关——这就是"最优子结构"。

本章用四个经典问题把这套方法练熟:钢条切割(一根长度为 \(n\) 的钢条怎么切卖得最多——可以类比为"一大笔订单怎么拆成几笔执行")、矩阵链乘法(一串矩阵相乘按什么顺序最省计算——直接用于因子风险模型)、最长公共子序列(两条序列最长的共同部分——可类比为比较两段价格形态)、最优二叉搜索树(把常查的东西放在离根近的位置)。每个问题都按"四步法"走:说清最优解长什么样 → 写出递推式 → 从小到大填表 → 从表里倒推出具体方案。

需要先想起来的数学。

  • 递推式(递归定义)。 用小规模的值定义大规模的值,比如 \(r_n=\max_i(p_i+r_{n-i})\),配上起点 \(r_0=0\)。你熟悉的例子是年金现值的递推:\(PV_n=\frac{C+PV_{n-1}}{1+r}\)。本章每个问题的核心都是写出这样一个式子。
  • \(\max\) 与 \(\min\) 记号。 \(\max_{1\le i\le n}f(i)\) 表示让 \(i\) 从 1 取到 \(n\),取 \(f(i)\) 的最大值;\(\arg\max\) 则表示取到最大值的那个 \(i\)。本章的 \(s\) 表、\(root\) 表存的就是 \(\arg\max\)/\(\arg\min\)。
  • 指数增长与多项式增长。 \(2^n\) 是指数:\(n\) 每加 1 翻一倍,\(n=40\) 时超过一万亿。\(n^2\)、\(n^3\) 是多项式:\(n=1000\) 时 \(n^3\) 是十亿,计算机几秒内可完成。动态规划的价值就是把前者变成后者。见 第 00 册第 04 章 级数与收敛。
  • 矩阵乘法的维数与计算量。 \(p\times q\) 矩阵乘 \(q\times r\) 矩阵得到 \(p\times r\) 矩阵;结果的每个元素是 \(q\) 个乘积之和,共 \(pr\) 个元素,所以要 \(pqr\) 次乘法。矩阵乘法满足结合律 \((AB)C=A(BC)\),但不满足交换律。见 第 00 册第 06 章 线性代数速成。
  • 反证法。 先假设结论不成立,推出矛盾。本章的"剪切—粘贴"论证全部是反证法:"如果子问题的解不是最优的,换成更优的,整体就更好,与整体最优矛盾。"见 第 00 册第 08 章 读懂数学证明与符号。

怎么读这一章。 15.1 节和 15.2 节(钢条切割)是全章最重要的部分,务必逐行读懂,包括"朴素递归为什么是指数时间"和两种实现。15.4 节"动态规划原理"是方法论的总结,必读,特别是 15.4.2 的"资源争用破坏最优子结构"——它直接关系到组合优化里的换手率约束能不能用 DP。15.3 节矩阵链乘法建议读懂递推式和原书例子,然后直接看 15.8.2 节的风险模型应用。15.5 节 LCS 读懂递推式即可。15.6 节最优 BST 第一次可以只看问题和递推式,算法细节可跳过。15.7 节思考题表格是下篇 15b 的地图,浏览即可。


15.1 动态规划是什么

动态规划和分治法一样,通过组合子问题的解来求解原问题(这里的 "programming" 指表格法,不是写程序)。区别在于:

  • 分治法把问题划分为互不相交的子问题,递归求解再合并;
  • 动态规划适用于子问题重叠的情形——不同的子问题共享子子问题。分治法会反复求解这些公共部分;动态规划对每个子子问题只求解一次,把答案存进表里,用到时查表。

白话解释:分治和动态规划的区别,用两个熟悉的例子对比。归并排序(分治)把 1000 个数分成左右两半各自排序,左半和右半没有任何共同元素,各算各的,不存在"重复劳动"。二叉树期权定价则不同:从根出发,"涨跌"和"跌涨"两条路走到同一个结点,"涨涨跌""涨跌涨""跌涨涨"三条路走到同一个结点。如果按路径逐条递归,同一个结点的价值会被算很多遍;动态规划的做法是把每个结点的价值算一次、记在表里,所有经过它的路径都来查表。"共享子子问题"指的就是这些被多条路径共用的结点。

动态规划通常用来求解最优化问题(optimization problem):问题可能有许多可行解,每个解有一个值,要找值最优的解。我们说求"一个最优解"而不是"唯一最优解",因为可能有多个解同时达到最优值。

设计动态规划算法的四个步骤:

  1. 刻画一个最优解的结构特征;
  2. 递归地定义最优解的值;
  3. 计算最优解的值,通常自底向上;
  4. 利用计算过程中记录的信息构造一个最优解。

前三步是基础。只要最优值、不要具体方案时可以省略第 4 步;要做第 4 步,常需在第 3 步顺带记录每个子问题上做的选择。


15.2 钢条切割

15.2.1 问题

Serling 公司把长钢条切成短钢条出售,切割不花钱。长度为 \(i\) 英寸的钢条售价 \(p_i\)。给定长度 \(n\) 和价格表,求使总收益 \(r_n\) 最大的切割方案。原书图 15.1 的价格表:

长度 \(i\) 1 2 3 4 5 6 7 8 9 10
价格 \(p_i\) 1 5 8 9 10 17 17 20 24 30

长度 \(n\) 的钢条在距左端 \(1,\dots,n-1\) 处都可以独立地选切或不切,共 \(2^{n-1}\) 种方案。\(n=4\) 时 8 种方案中切成 \(2+2\) 收益最大,为 10。前 10 个长度的最优收益是:\(r_1=1\),\(r_2=5\),\(r_3=8\),\(r_4=10\,(2+2)\),\(r_5=13\,(2+3)\),\(r_6=17\)(不切),\(r_7=18\,(1+6\) 或 \(2+2+3)\),\(r_8=22\,(2+6)\),\(r_9=25\,(3+6)\),\(r_{10}=30\)(不切)。

15.2.2 最优子结构与递归式

第一种看法:要么不切,要么先切成长 \(i\) 和 \(n-i\) 两段,再对两段分别最优切割:

\[r_n=\max(p_n,\ r_1+r_{n-1},\ r_2+r_{n-2},\ \dots,\ r_{n-1}+r_1).\]
原问题的最优解由两个子问题的最优解组成,而且两个子问题可以独立求解——称这个问题具有最优子结构(optimal substructure)。

更简单的看法:把方案看作"从左端切下第一段长度 \(i\)(这一段不再切),再对剩下的 \(n-i\) 继续切"。不切就是 \(i=n\)、剩余 0、\(r_0=0\):

\[r_n=\max_{1\le i\le n}\,(p_i+r_{n-i}).\tag{15.2}\]
这样最优解只涉及一个子问题(剩余部分)。

推导拆解:用价格表手算 \(r_4\),看式 (15.2) 怎么用。第一段切多长有 4 种选择: 切 1 段长 1:\(p_1+r_3=1+8=9\); 切 1 段长 2:\(p_2+r_2=5+5=10\); 切 1 段长 3:\(p_3+r_1=8+1=9\); 整根不切:\(p_4+r_0=9+0=9\)。 取最大,\(r_4=10\)。关键在于:算 \(r_4\) 时用到的 \(r_3,r_2,r_1\) 都是已经算好的最优值,我们不关心长度 3 的那段内部怎么切,只要知道它最多能卖 8。这就是"最优子结构"的含义:剩下那段无论怎么来的,都应该按它自己的最优方式切。 交易类比:一笔 4 万股的大单,第一笔先下多少(1、2、3 还是 4 万股)?下完第一笔后,剩下的数量也应该按"剩余数量的最优拆法"执行。15b 章的拆单模型就是这个结构,只是"价格表"换成了冲击成本函数。

15.2.3 朴素递归为什么是指数时间

def CUT_ROD(p, n):
    if n == 0: return 0
    q = -inf
    for i in range(1, n + 1):
        q = max(q, p[i] + CUT_ROD(p, n - i))
    return q

它对 \(j=0,1,\dots,n-1\) 的每一个都发起递归调用,而这些调用又各自重复同样的事。设 \(T(n)\) 为以 \(n\) 为参数时的调用总数(含自身),则

\[T(0)=1,\qquad T(n)=1+\sum_{j=0}^{n-1}T(j).\]
由 \(T(n)-T(n-1)=T(n-1)\) 得 \(T(n)=2^n\)。事后看并不奇怪:它显式考察了全部 \(2^{n-1}\) 种切法,递归树的 \(2^{n-1}\) 个叶对应 \(2^{n-1}\) 种方案。原书说 \(n=40\) 时很可能要跑一个多小时,\(n\) 每加 1 时间约翻一倍。

推导拆解:\(T(n)=2^n\) 的两行推导。写出 \(T(n)=1+T(0)+T(1)+\cdots+T(n-1)\) 和 \(T(n-1)=1+T(0)+\cdots+T(n-2)\),两式相减,右边只剩 \(T(n-1)\),所以 \(T(n)-T(n-1)=T(n-1)\),即 \(T(n)=2T(n-1)\)。从 \(T(0)=1\) 起每步翻倍,\(T(n)=2^n\)。 浪费在哪里:算 CUT_ROD(4) 时要调用 CUT_ROD(3)、(2)、(1)、(0);而 CUT_ROD(3) 内部又调用 (2)、(1)、(0);CUT_ROD(2) 内部又调用 (1)、(0)……CUT_ROD(1) 被调用了 4 次,CUT_ROD(0) 被调用了 8 次,每次都从头重算一遍同样的答案。这就像每次要用"5 年期年金现值系数"时都从头逐年贴现一遍,而不是查一张已经算好的年金系数表。

15.2.4 两种动态规划实现

问题出在重复求解相同的子问题。动态规划的办法是:每个子问题只解一次并存下结果。这是一种时空权衡(time-memory trade-off):用额外的内存换时间,能把指数时间降为多项式时间——前提是不同子问题的个数是多项式的,且每个子问题能在多项式时间内求解。

带备忘的自顶向下法(top-down with memoization):保持递归的写法,但把每个子问题的结果存在数组或散列表里;进入过程先查表,查到就直接返回。(memoization 来自 memo,不是 memorization 的笔误。)

def MEMOIZED_CUT_ROD(p, n):
    r = [-inf] * (n + 1)              # -inf 表示"尚未计算"(收益总是非负)
    return MEMOIZED_CUT_ROD_AUX(p, n, r)

def MEMOIZED_CUT_ROD_AUX(p, n, r):
    if r[n] >= 0: return r[n]
    if n == 0: q = 0
    else:
        q = -inf
        for i in range(1, n + 1):
            q = max(q, p[i] + MEMOIZED_CUT_ROD_AUX(p, n - i, r))
    r[n] = q
    return q

自底向上法(bottom-up method):定义子问题的"规模",使每个子问题只依赖更小的子问题;按规模从小到大求解,求解某个子问题时它依赖的子问题都已解好。

def BOTTOM_UP_CUT_ROD(p, n):
    r = [0] * (n + 1)
    for j in range(1, n + 1):
        q = -inf
        for i in range(1, j + 1):
            q = max(q, p[i] + r[j - i])
        r[j] = q
    return r[n]

两者都是 \(\Theta(n^2)\) 时间:自底向上的双重循环迭代次数构成等差数列;备忘法中每个规模的子问题只真正求解一次,求解规模 \(j\) 时 for 循环迭代 \(j\) 次,总数同样是等差数列。自底向上没有递归调用开销,常数因子通常更小;备忘法的优势是只求解真正用到的子问题。

白话解释:两种实现的区别是"按需记账"还是"提前造表"。 备忘法像一个做估值的分析师:老板问"10 年期怎么切最优",他想"那我得先知道 9 年、8 年……的答案",于是一层层往下问;每算出一个答案就记到小本子上,下次有人问同一个问题直接翻本子。代码里的 r[n] >= 0 就是"翻本子看有没有记过"。 自底向上像编制年金系数表:先填 1 年期,再填 2 年期(用到 1 年期),再填 3 年期(用到 1、2 年期)……填到 10 年期时,所需的一切都已在表上。代码里外层循环 for j in 1..n 就是"按期限从短到长逐行填表"。 两者算的东西完全一样,数字也一样,只是顺序不同。为什么都是 \(\Theta(n^2)\):长度 \(j\) 的子问题要试 \(j\) 种第一刀,\(1+2+\cdots+n=\frac{n(n+1)}2\)。

15.2.5 子问题图

**子问题图(subproblem graph)**是一个有向图:每个不同的子问题一个顶点,若求解 \(x\) 时要直接用到 \(y\) 的最优解,就画一条边 \((x,y)\)。它是递归树的"压缩版"——把递归树中代表同一子问题的结点合并。钢条切割 \(n=4\) 的子问题图有顶点 \(0,\dots,4\),每个顶点指向所有更小的顶点(原书图 15.4)。

  • 自底向上法按子问题图的逆拓扑序处理顶点:一个子问题在它依赖的子问题都解完之前不会被考虑。
  • 带备忘的自顶向下法相当于在子问题图上做深度优先搜索。
  • 运行时间通常与子问题图的顶点数加边数成线性:每个子问题只解一次,解它的时间通常正比于它的出度。

15.2.6 重构解

只返回最优值还不够,还要知道怎么切。对每个子问题额外记录"导致最优值的选择":

def EXTENDED_BOTTOM_UP_CUT_ROD(p, n):
    r, s = [0] * (n + 1), [0] * (n + 1)   # s[j]:长度 j 的最优方案中第一段的长度
    for j in range(1, n + 1):
        q = -inf
        for i in range(1, j + 1):
            if q < p[i] + r[j - i]:
                q, s[j] = p[i] + r[j - i], i
        r[j] = q
    return r, s

然后反复输出 \(s[n]\) 并令 \(n\leftarrow n-s[n]\)。原书样例的结果:

\(i\) 0 1 2 3 4 5 6 7 8 9 10
\(r[i]\) 0 1 5 8 10 13 17 18 22 25 30
\(s[i]\) 0 1 2 3 2 2 6 1 2 3 10

\(n=10\) 时只输出 10(不切);\(n=7\) 时输出 1 和 6。

两个变体值得注意:习题 15.1-2 说明"按单位长度价格 \(p_i/i\) 最高的贪心"不一定最优(例如 \(p_1=1,p_2=5,p_3=8\)、\(n=4\):贪心先切密度最大的长度 3,得 \(8+1=9<10\));习题 15.1-3 给每次切割加固定成本 \(c\),递归式变为 \(r_n=\max\big(p_n,\ \max_{1\le i<n}(p_i+r_{n-i}-c)\big)\)——这和"每次调仓付固定费用"的决策结构一模一样。


15.3 矩阵链乘法

15.3.1 问题

给定矩阵链 \(\langle A_1,A_2,\dots,A_n\rangle\),\(A_i\) 的维数是 \(p_{i-1}\times p_i\),求乘积 \(A_1A_2\cdots A_n\)。矩阵乘法满足结合律,所有加括号方式结果相同,但代价可能相差悬殊。一个 \(p\times q\) 矩阵乘 \(q\times r\) 矩阵需要 \(pqr\) 次标量乘法。

原书例:\(A_1:10\times100\),\(A_2:100\times5\),\(A_3:5\times50\)。按 \(((A_1A_2)A_3)\) 共 \(5000+2500=7500\) 次;按 \((A_1(A_2A_3))\) 共 \(25000+50000=75000\) 次,慢 10 倍。

矩阵链乘法问题:求使标量乘法次数最少的完全括号化方案。注意目标不是真的去乘,而是确定最便宜的计算顺序。

方案数 \(P(n)\) 满足 \(P(1)=1\),\(P(n)=\sum_{k=1}^{n-1}P(k)P(n-k)\),其解是 Catalan 数,增长为 \(\Omega(4^n/n^{3/2})\)。穷举不可行。

15.3.2 四步法

步骤 1:最优括号化的结构。 记 \(A_{i..j}=A_iA_{i+1}\cdots A_j\)。任何 \(i<j\) 的括号化都要在某个 \(k\)(\(i\le k<j\))处把链分开:先算 \(A_{i..k}\) 和 \(A_{k+1..j}\),再相乘。最优子结构:若最优方案在 \(k\) 处分开,则前缀 \(A_i\cdots A_k\) 的括号化必须是该子链的最优括号化——否则换成更优的方案,就得到更便宜的整体方案,矛盾。后缀同理。

步骤 2:递归解。 令 \(m[i,j]\) 为计算 \(A_{i..j}\) 的最少标量乘法次数:

\[m[i,j]=\begin{cases}0 & i=j,\\ \min_{i\le k<j}\{m[i,k]+m[k+1,j]+p_{i-1}p_kp_j\} & i<j.\end{cases}\tag{15.7}\]
记 \(s[i,j]\) 为取得最优值的分割点 \(k\)。

推导拆解:式 (15.7) 的三项分别是什么。在 \(k\) 处分开时,计算 \(A_{i..j}\) 分三件事:先算左半 \(A_{i..k}\),最少要 \(m[i,k]\) 次;再算右半 \(A_{k+1..j}\),最少要 \(m[k+1,j]\) 次;最后把两个结果相乘。左半结果的维数是"\(A_i\) 的行数 × \(A_k\) 的列数" \(=p_{i-1}\times p_k\),右半结果是 \(p_k\times p_j\),相乘要 \(p_{i-1}p_kp_j\) 次。对所有可能的 \(k\) 取最小即可。 用原书三矩阵例子核对:\(p=\langle10,100,5,50\rangle\)。\(k=1\):\(m[1,1]+m[2,3]+p_0p_1p_3=0+25000+10\cdot100\cdot50=75000\);\(k=2\):\(m[1,2]+m[3,3]+p_0p_2p_3=5000+0+10\cdot5\cdot50=7500\)。取 \(k=2\),即 \(((A_1A_2)A_3)\)。 为什么子问题要用两端都可变的 \(i..j\):在 \(k\) 处一刀切下去,产生的是"前段"和"后段",后段 \(A_{k+1..j}\) 不是从第 1 个矩阵开始的,所以单用"前缀"描述不了它。

步骤 3:计算最优代价。 不同子问题只有 \(\binom n2+n=\Theta(n^2)\) 个,但直接递归会在递归树的不同分支里反复遇到同一子问题——这就是重叠子问题(overlapping subproblems)。由于 \(m[i,j]\) 只依赖更短的链,按链长 \(l=j-i+1\) 递增的顺序填表:

def MATRIX_CHAIN_ORDER(p):
    n = len(p) - 1
    m = [[0] * (n + 1) for _ in range(n + 1)]
    s = [[0] * (n + 1) for _ in range(n + 1)]
    for l in range(2, n + 1):                 # 链长
        for i in range(1, n - l + 2):
            j = i + l - 1
            m[i][j] = inf
            for k in range(i, j):
                q = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j]
                if q < m[i][j]:
                    m[i][j], s[i][j] = q, k
    return m, s

原书图 15.5 例:\(p=\langle30,35,15,5,10,20,25\rangle\)。表中各值为 \(m[1,2]=15750\),\(m[2,3]=2625\),\(m[3,4]=750\),\(m[4,5]=1000\),\(m[5,6]=5000\);\(m[1,3]=7875\),\(m[2,4]=4375\),\(m[3,5]=2500\),\(m[4,6]=3500\);\(m[1,4]=9375\),\(m[2,5]=7125\),\(m[3,6]=5375\);\(m[1,5]=11875\),\(m[2,6]=10500\);\(m[1,6]=15125\)。以 \(m[2,5]\) 为例:

\[m[2,5]=\min\begin{cases}m[2,2]+m[3,5]+p_1p_2p_5=0+2500+35\cdot15\cdot20=13000,\\ m[2,3]+m[4,5]+p_1p_3p_5=2625+1000+35\cdot5\cdot20=7125,\\ m[2,4]+m[5,5]+p_1p_4p_5=4375+0+35\cdot10\cdot20=11375\end{cases}=7125.\]

三重循环,时间 \(\Theta(n^3)\)(习题 15.2-5 证明了下界),空间 \(\Theta(n^2)\)。

步骤 4:构造最优解。 \(s[1,n]\) 告诉我们最后一次乘法在哪里分开,再递归处理两边:

def PRINT_OPTIMAL_PARENS(s, i, j):
    if i == j: print(f"A{i}", end="")
    else:
        print("(", end="")
        PRINT_OPTIMAL_PARENS(s, i, s[i][j])
        PRINT_OPTIMAL_PARENS(s, s[i][j] + 1, j)
        print(")", end="")

图 15.5 的例子输出 \(((A_1(A_2A_3))((A_4A_5)A_6))\)。


15.4 动态规划原理

从工程角度看,什么时候该想到动态规划?两个要素:最优子结构和重叠子问题。

15.4.1 最优子结构

如果问题的最优解包含其子问题的最优解,就说它具有最优子结构。发掘最优子结构有一个通用模式:

  1. 证明问题的解要做一个选择(钢条的第一刀切在哪里、矩阵链在哪里分开),做出选择后留下一个或多个子问题;
  2. 假定已经知道哪个选择会导向最优解——先不管怎么知道的;
  3. 给定这个选择,确定会产生哪些子问题,以及如何最好地刻画子问题空间;
  4. 用**剪切—粘贴(cut-and-paste)**论证子问题的解必须是最优的:如果某个子问题的解不是最优的,把它"剪掉"、"粘上"该子问题的最优解,就得到原问题更好的解,与原解最优矛盾。

金融直觉:剪切—粘贴论证在金融里叫贝尔曼最优性原理,你在美式期权倒推定价中已经默默用过:如果从 0 时刻看的最优行权策略在某个中间结点上"继续持有",那么从这个结点往后看,这个策略也必须是该结点上的最优策略。否则把这个结点之后的部分"剪掉",换成该结点上真正最优的策略"粘上",整体价值就更高,与"原策略最优"矛盾。这句话成立的前提是:结点之后的最优决策只取决于当前状态(股价、剩余期限),不取决于你是如何走到这里的。一旦未来决策依赖于过去(例如路径依赖的亚式期权要看历史均价),就得把相关的历史信息(如累计均价)放进状态里,否则最优子结构不成立——这正是 15.4.2 节要讲的问题。

刻画子问题空间的经验:尽量简单,必要时再扩展。钢条切割用"长度为 \(i\) 的钢条"就够了。矩阵链若只用前缀 \(A_1\cdots A_j\) 作子问题是不够的——在 \(k\) 处分开会产生 \(A_{k+1}\cdots A_j\),所以子问题必须两端都可变。

两类差异决定了不同问题的 DP 形态:最优解用到几个子问题(钢条 1 个,矩阵链 2 个);在多少种选择中挑选(钢条 \(n\) 种,矩阵链 \(j-i\) 种)。运行时间 ≈ 子问题总数 × 每个子问题的选择数:钢条 \(\Theta(n)\times O(n)=O(n^2)\),矩阵链 \(\Theta(n^2)\times O(n)=O(n^3)\)。用子问题图同样可以得到这个估计。原解的代价通常是子问题代价加上选择本身直接带来的代价(钢条中的 \(p_i\),矩阵链中的 \(p_{i-1}p_kp_j\))。

与贪心算法的区别(下一章):贪心问题也有最优子结构,但贪心不先求子问题再选择,而是先做"当前看来最好"的选择,再解剩下的唯一子问题。

15.4.2 最优子结构不总成立:最短路径与最长简单路径

在有向图中:

  • 无权最短路径:从 \(u\) 到 \(v\) 边数最少的路径。它有最优子结构:若 \(w\) 在最短路径 \(p\) 上,把 \(p\) 分成 \(u\leadsto w\) 与 \(w\leadsto v\),两段都必须是最短的(剪切—粘贴)。
  • 无权最长简单路径:边数最多的简单路径。它没有最优子结构。原书图 15.6:顶点 \(q,r,s,t\),\(q\to r\to t\) 是 \(q\) 到 \(t\) 的最长简单路径,但 \(q\to r\) 不是 \(q\) 到 \(r\) 的最长简单路径(\(q\to s\to t\to r\) 更长);而把两个子问题的最长路径拼起来,\(q\to s\to t\to r\to q\to s\to t\) 甚至不是简单路径。这个问题是 NP 完全的。

根本原因是子问题是否独立(independent):一个子问题的解是否影响另一个子问题的解。最长简单路径的两个子问题争用顶点——第一个子问题用掉了 \(s\) 和 \(t\),第二个就不能再用。最短路径的子问题天然不共享顶点(可以证明:若两段最短路径除 \(w\) 外还共享顶点 \(x\),就能剪掉 \(x\to w\to x\) 这一圈,得到更短的路径,矛盾)。钢条切割和矩阵链的子问题也都独立。

这个"资源争用破坏最优子结构"的观点在量化里很重要。原书习题 15.3-5 说:若每种长度的钢条数量有上限 \(l_i\),最优子结构就不成立——两个子问题共享"配额"。同样,思考题 15-10(d) 说:若规定任何时刻单一投资不得超过某个金额,多期投资问题就失去最优子结构。组合优化中的总量约束、换手率约束都是这样的共享资源,往往需要把资源本身放进状态(状态维数上升),或改用数学规划。

15.4.3 重叠子问题

第二个要素是子问题空间要"足够小":递归算法反复求解相同的子问题,而不是不断产生新问题。适合分治法的问题通常每一步都产生全新的子问题(比如归并排序),备忘对它们没有帮助(习题 15.3-2)。

白话解释:判断一个问题能否用 DP,可以问两个问题。第一,"做完第一个决定后,剩下的问题是不是和原问题同一类、只是规模更小?"(最优子结构)钢条切掉第一段后剩下的还是一根钢条;拆单下完第一笔后剩下的还是一笔待执行的单。第二,"不同的决策路线会不会走到同一个剩余问题?"(重叠子问题)先切 1 再切 2、先切 2 再切 1,剩下的都是长度 \(n-3\) 的钢条;先卖 1 万再卖 2 万、先卖 2 万再卖 1 万,剩下的都是"还有 \(Q-3\) 万股"。两个答案都是"是",DP 就有用武之地。归并排序只满足第一条,左右两半永不相同,记下答案也不会被再用到,所以它只是分治。

"独立"与"重叠"听起来矛盾,其实是两个不同维度:同一问题的两个子问题不共享资源叫独立;同一个子问题作为不同问题的子问题反复出现叫重叠。

看直接按式 15.7 写的递归:

def RECURSIVE_MATRIX_CHAIN(p, i, j):
    if i == j: return 0
    m[i][j] = inf
    for k in range(i, j):
        q = (RECURSIVE_MATRIX_CHAIN(p, i, k) + RECURSIVE_MATRIX_CHAIN(p, k + 1, j)
             + p[i - 1] * p[k] * p[j])
        m[i][j] = min(m[i][j], q)
    return m[i][j]

设 \(T(n)\) 为其时间,则 \(T(1)\ge1\),\(T(n)\ge1+\sum_{k=1}^{n-1}(T(k)+T(n-k)+1)\),即

\[T(n)\ge2\sum_{i=1}^{n-1}T(i)+n.\]
用代入法证 \(T(n)\ge2^{n-1}\):\(T(n)\ge2\sum_{i=1}^{n-1}2^{i-1}+n=2(2^{n-1}-1)+n=2^n-2+n\ge2^{n-1}\)。所以它至少是指数时间,而不同子问题只有 \(\Theta(n^2)\) 个。

15.4.4 重构解与备忘

实践中总是把每个子问题的选择存进表(如 \(s[i,j]\)),否则从代价表里反推选择要 \(\Theta(j-i)\) 时间。

给上面的递归加备忘,就得到 \(O(n^3)\) 的 MEMOIZED-MATRIX-CHAIN:表项初值设为 \(\infty\) 表示"未计算",第一次遇到时计算并存表,之后直接返回。复杂度分析:表项共 \(\Theta(n^2)\) 个,每个只"真正计算"一次,每次计算发出 \(O(n)\) 个递归调用;直接返回的调用共 \(O(n^3)\) 次、每次 \(O(1)\)。总计 \(O(n^3)\)。

怎么选? 若所有子问题都至少要解一次,自底向上通常快一个常数因子,而且便于利用表访问的规律压缩空间(如 LCS 只保留两行)。若有些子问题根本不需要解,备忘法只解需要的,更有优势。另外,当子问题的参数不是简单下标时(例如状态是一个元组),备忘法用散列表存储最自然。


15.5 最长公共子序列

15.5.1 问题

序列 \(Z=\langle z_1,\dots,z_k\rangle\) 是 \(X=\langle x_1,\dots,x_m\rangle\) 的子序列(subsequence),若存在严格递增的下标 \(i_1<\dots<i_k\) 使 \(x_{i_j}=z_j\)(即删去若干元素,不必连续)。\(Z\) 同时是 \(X\) 和 \(Y\) 的子序列时称为公共子序列,最长的称为最长公共子序列(LCS)。

原书例:\(X=\langle A,B,C,B,D,A,B\rangle\),\(Y=\langle B,D,C,A,B,A\rangle\)。\(\langle B,C,A\rangle\) 是公共子序列但不是最长的;\(\langle B,C,B,A\rangle\) 和 \(\langle B,D,A,B\rangle\) 都是 LCS,长度 4。

暴力枚举 \(X\) 的 \(2^m\) 个子序列是指数时间。记前缀 \(X_i=\langle x_1,\dots,x_i\rangle\),子问题自然对应前缀对。

15.5.2 最优子结构

定理 15.1:设 \(Z=\langle z_1,\dots,z_k\rangle\) 是 \(X\) 与 \(Y\) 的任一 LCS。

  1. 若 \(x_m=y_n\),则 \(z_k=x_m=y_n\),且 \(Z_{k-1}\) 是 \(X_{m-1}\) 与 \(Y_{n-1}\) 的一个 LCS;
  2. 若 \(x_m\ne y_n\),则 \(z_k\ne x_m\) 意味着 \(Z\) 是 \(X_{m-1}\) 与 \(Y\) 的一个 LCS;
  3. 若 \(x_m\ne y_n\),则 \(z_k\ne y_n\) 意味着 \(Z\) 是 \(X\) 与 \(Y_{n-1}\) 的一个 LCS。

证明. (1) 若 \(z_k\ne x_m\),把 \(x_m=y_n\) 追加到 \(Z\) 后得到更长的公共子序列,矛盾,所以 \(z_k=x_m\)。\(Z_{k-1}\) 是 \(X_{m-1},Y_{n-1}\) 的长 \(k-1\) 的公共子序列;若有更长的 \(W\),追加 \(x_m\) 就得到 \(X,Y\) 长度超过 \(k\) 的公共子序列,矛盾。(2) \(z_k\ne x_m\) 时 \(Z\) 是 \(X_{m-1}\) 与 \(Y\) 的公共子序列;若有更长的 \(W\),它也是 \(X\) 与 \(Y\) 的公共子序列,与 \(Z\) 最长矛盾。(3) 对称。\(\square\)

15.5.3 递归式与算法

令 \(c[i,j]\) 为 \(X_i\) 与 \(Y_j\) 的 LCS 长度:

\[c[i,j]=\begin{cases}0 & i=0\ \text{或}\ j=0,\\ c[i-1,j-1]+1 & i,j>0,\ x_i=y_j,\\ \max(c[i,j-1],\,c[i-1,j]) & i,j>0,\ x_i\ne y_j.\end{cases}\tag{15.9}\]

白话解释:式 (15.9) 只做一件事:比较两条序列的最后一个元素。 如果两个最后元素相同(比如都是 B),那它一定可以作为公共子序列的最后一个,答案 = 两条序列各自去掉末尾后的 LCS 长度 + 1。 如果不同,那这两个末尾元素不可能同时用上,至少有一个是多余的:要么扔掉 \(X\) 的末尾,要么扔掉 \(Y\) 的末尾,两种情况各算一次取较大者。 每一步都让至少一条序列变短,最后变成空序列,答案为 0。子问题就是"\(X\) 的前 \(i\) 个与 \(Y\) 的前 \(j\) 个",共 \((m+1)(n+1)\) 个,填成一张表:行是 \(i\),列是 \(j\),每格只看它左边、上边、左上三个格子。 量化里的类比:把两只股票的日涨跌编码成"涨/跌/平"序列,LCS 长度衡量两者走势的"共同骨架",允许中间有若干天对不上。

注意问题条件会限制要考虑的子问题:\(x_i=y_j\) 时只看 \(c[i-1,j-1]\),否则看另两个。共 \(\Theta(mn)\) 个子问题,按行主序填表,每项 \(O(1)\):

def LCS_LENGTH(X, Y):
    m, n = len(X), len(Y)
    c = [[0] * (n + 1) for _ in range(m + 1)]
    b = [[None] * (n + 1) for _ in range(m + 1)]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if X[i] == Y[j]:
                c[i][j] = c[i - 1][j - 1] + 1; b[i][j] = "↖"
            elif c[i - 1][j] >= c[i][j - 1]:
                c[i][j] = c[i - 1][j];         b[i][j] = "↑"
            else:
                c[i][j] = c[i][j - 1];         b[i][j] = "←"
    return c, b

时间和空间都是 \(\Theta(mn)\)。从 \(b[m,n]\) 沿箭头回溯,遇到"↖"就输出 \(x_i\),递归地按正序打印,\(O(m+n)\)。原书图 15.8 输出 BCBA。

改进:\(b\) 表可以去掉,因为 \(c[i,j]\) 只依赖三个相邻表项,可 \(O(1)\) 判断当初用了哪个(习题 15.4-2);若只要长度,只需保留两行,空间降到 \(O(\min(m,n))\)(习题 15.4-4),但这样就无法 \(O(m+n)\) 回溯出 LCS 本身。

相关问题:最长单调递增子序列(LIS),\(O(n^2)\) 可以化为原序列与其排序结果的 LCS;用"各长度候选子序列的最小末元素"加二分查找可做到 \(O(n\lg n)\)(习题 15.4-5、15.4-6)。编辑距离(思考题 15-5)与 LCS 同构,只是把"匹配加 1"换成各种操作的代价;量化里常用的动态时间规整(DTW)——比较两段节奏不同的价格形态——也是同一类二维表 DP。


15.6 最优二叉搜索树

15.6.1 问题

给定 \(n\) 个有序关键字 \(k_1<\dots<k_n\),搜索 \(k_i\) 的概率为 \(p_i\);还有 \(n+1\) 个伪关键字 \(d_0,\dots,d_n\) 表示落在关键字之间(及两端之外)的不成功搜索,概率为 \(q_i\),\(\sum p_i+\sum q_i=1\)。搜索代价为访问的结点数,即深度加 1。树 \(T\) 的期望搜索代价为

\[E[\text{代价}]=1+\sum_{i=1}^n\text{depth}_T(k_i)\,p_i+\sum_{i=0}^n\text{depth}_T(d_i)\,q_i.\tag{15.11}\]
目标是构造期望代价最小的 BST。动机是英译法词典:常用词应靠近根。

原书图 15.9 例(\(n=5\)):\(p=(0.15,0.10,0.05,0.10,0.20)\),\(q=(0.05,0.10,0.05,0.05,0.05,0.10)\)。树 (a)(根 \(k_2\))期望代价 2.80;树 (b) 期望代价 2.75,是最优的。注意两点:最优树不一定高度最小;也不能把概率最大的关键字放在根——\(k_5\) 概率最大,但以 \(k_5\) 为根的最好的树期望代价是 2.85。

15.6.2 结构与递归式

BST 的任一子树都包含连续的关键字 \(k_i,\dots,k_j\),叶为 \(d_{i-1},\dots,d_j\)。若最优树有包含 \(k_i..k_j\) 的子树,它对这个子问题也必须最优(剪切—粘贴)。在 \(k_i..k_j\) 中选一个根 \(k_r\),左子树为 \(k_i..k_{r-1}\) 的最优树,右子树为 \(k_{r+1}..k_j\) 的最优树。

一棵子树挂到某个结点下时,其中每个结点深度加 1,期望代价增加子树的总概率

\[w(i,j)=\sum_{l=i}^jp_l+\sum_{l=i-1}^jq_l.\]
于是以 \(k_r\) 为根时 \(e[i,j]=p_r+(e[i,r-1]+w(i,r-1))+(e[r+1,j]+w(r+1,j))=e[i,r-1]+e[r+1,j]+w(i,j)\),
\[e[i,j]=\begin{cases}q_{i-1} & j=i-1,\\ \min_{i\le r\le j}\{e[i,r-1]+e[r+1,j]+w(i,j)\} & i\le j.\end{cases}\tag{15.14}\]

推导拆解:为什么挂到下一层就要"加上整棵子树的总概率"?搜索代价 = 访问的结点数。一棵子树单独看时,它的根在第 0 层;挂到 \(k_r\) 下面之后,它里面每个结点都往下移了一层,每次搜索到其中任何结点都要多访问 1 个结点(就是 \(k_r\))。期望代价增加 \(1\times\)(落在这棵子树的概率之和)\(=w\)。 合并的那一步:\(p_r+w(i,r-1)+w(r+1,j)\) 恰好是"根 \(k_r\) 的概率 + 左子树总概率 + 右子树总概率",也就是 \(k_i..k_j\) 这一段(含伪关键字)的总概率 \(w(i,j)\)。所以无论选哪个根,加上的那一项都是同一个 \(w(i,j)\),只有 \(e[i,r-1]+e[r+1,j]\) 随 \(r\) 变化。 交易类比:行情系统把最常被查询的代码放在离"入口"最近的位置;但最优结构不是简单的"最热门的放根上",因为根的选择还决定了左右两边如何分割,要在全局上权衡。

15.6.3 算法

用 \(w[i,j]=w[i,j-1]+p_j+q_j\) 递推 \(w\),按子树关键字个数 \(l\) 递增填表,最内层尝试每个候选根:

def OPTIMAL_BST(p, q, n):
    e = [[0] * (n + 1) for _ in range(n + 2)]       # e[1..n+1][0..n]
    w = [[0] * (n + 1) for _ in range(n + 2)]
    root = [[0] * (n + 1) for _ in range(n + 1)]
    for i in range(1, n + 2):
        e[i][i - 1] = w[i][i - 1] = q[i - 1]
    for l in range(1, n + 1):
        for i in range(1, n - l + 2):
            j = i + l - 1
            e[i][j] = inf
            w[i][j] = w[i][j - 1] + p[j] + q[j]
            for r in range(i, j + 1):
                t = e[i][r - 1] + e[r + 1][j] + w[i][j]
                if t < e[i][j]:
                    e[i][j], root[i][j] = t, r
    return e, root

时间 \(\Theta(n^3)\),空间 \(\Theta(n^2)\)。对图 15.9 的分布,\(e[1,5]=2.75\),\(w[1,5]=1.00\);根表中 \(root[1,5]=2\)、\(root[2,5]=4\)、\(root[3,5]=5\)。最优树的结构(习题 15.5-1):\(k_2\) 为根,左孩子 \(k_1\)(其下为 \(d_0,d_1\)),右孩子 \(k_5\);\(k_5\) 的左孩子 \(k_4\)、右孩子 \(d_5\);\(k_4\) 的左孩子 \(k_3\)(其下为 \(d_2,d_3\))、右孩子 \(d_4\)。

(勘误说明:本册所依据的精读笔记把根表的一项记作 "\(root[2,5]=5\)";按递归式实际计算,\(root[2,5]=4\),等于 5 的是 \(root[3,5]\)——这也与上面的树结构一致:包含 \(k_3..k_5\) 的右子树以 \(k_5\) 为根。下文代码给出了全部数值。)

Knuth 优化(习题 15.5-4):总存在最优根满足 \(root[i,j-1]\le root[i,j]\le root[i+1,j]\)。利用这个单调性,最内层循环只需在一个窗口内搜索,每条对角线上的总工作量为 \(O(n)\),整体降到 \(\Theta(n^2)\)。这类"决策单调性"技巧在最优分段等问题中也常用。


15.7 原书思考题速览

原书第 15 章的 12 道思考题几乎每一道都对应一个经典 DP 模型,其中好几道与量化直接相关:

题号 问题 DP 要点 量化对应
15-1 DAG 上的最长带权路径 按拓扑序 DP,\(O(V+E)\);无环所以没有资源争用 事件依赖图上的关键路径
15-2 最长回文子序列 区间 DP,\(O(n^2)\) —
15-3 双调欧几里得旅行商 从左到右扫描,\(O(n^2)\) —
15-4 整齐打印 分段代价为行末空格立方和,\(O(n^2)\) 时间序列最优分段/变点检测的同构问题
15-5 编辑距离与序列比对 二维表,\(\Theta(mn)\) 形态匹配、DTW
15-6 公司聚会 树形 DP,每个结点算"来/不来"两个值,\(O(n)\) —
15-7 Viterbi 算法 按时间逐层取最可能前驱,\(O(k(V+E))\) HMM 市场状态解码(见 15b)
15-8 接缝裁剪 \(D[i,j]=d[i,j]+\min(D[i-1,j-1],D[i-1,j],D[i-1,j+1])\) —
15-9 拆分字符串 区间 DP,类似矩阵链 —
15-10 投资策略规划 状态"第 \(j\) 年持有投资 \(i\)";加持仓上限后失去最优子结构 带换仓费用的多期轮动(见 15b)
15-11 库存规划 状态"月份 × 月末库存" 拆单执行的同构问题(见 15b)
15-12 签约自由球员 分组背包,预算为状态 预算约束下的策略/因子选择

习题 15.3-6 的货币兑换也值得一提:\(n\) 种货币按汇率 \(r_{ij}\) 兑换,佣金为零时,最优兑换序列具有最优子结构;佣金随兑换次数变化时则不一定。把汇率取负对数后,"乘积最大"变成"路径和最小",寻找套利环就是在图上找负权环——见本册第 24 章的 Bellman-Ford 算法。


15.8 量化实战

15.8.1 验证四个经典问题

先用代码把本章的四个经典例子全部跑一遍,同时对比钢条切割三种实现的速度:

import numpy as np, time
INF = float("inf")

# ---------- 钢条切割:朴素递归 vs 备忘 vs 自底向上 ----------
p = [0, 1, 5, 8, 9, 10, 17, 17, 20, 24, 30]            # 原书图 15.1 价格表
def cut_rod(p, n):                                      # 朴素递归,2^n 次调用
    if n == 0: return 0
    return max(p[i] + cut_rod(p, n - i) for i in range(1, n + 1))
def memoized_cut_rod(p, n):                             # 带备忘的自顶向下:用散列表(dict)存子问题解
    memo = {0: 0}
    def aux(j):
        if j not in memo:
            memo[j] = max(p[i] + aux(j - i) for i in range(1, j + 1))
        return memo[j]
    return aux(n)
def extended_bottom_up_cut_rod(p, n):
    r, s = [0] * (n + 1), [0] * (n + 1)
    for j in range(1, n + 1):
        q = -INF
        for i in range(1, j + 1):
            if q < p[i] + r[j - i]: q, s[j] = p[i] + r[j - i], i
        r[j] = q
    return r, s
r, s = extended_bottom_up_cut_rod(p, 10)
print("r =", r); print("s =", s)
n, cuts = 7, []
while n > 0: cuts.append(s[n]); n -= s[n]
print("长度 7 的最优切法:", cuts)
pp = p + [p[10] + 0] * 30                               # 把价格表延长,用来计时
for n in [16, 20, 22]:
    t0 = time.perf_counter(); a = cut_rod(pp, n); t1 = time.perf_counter()
    b = memoized_cut_rod(pp, n); t2 = time.perf_counter()
    c = extended_bottom_up_cut_rod(pp, n)[0][n]; t3 = time.perf_counter()
    print(f"n={n}: 朴素递归 {t1-t0:6.3f}s | 备忘 {1e3*(t2-t1):6.3f}ms | 自底向上 {1e3*(t3-t2):6.3f}ms | 三者相等 {a == b == c}")

# ---------- 矩阵链乘法:原书图 15.5 ----------
def matrix_chain_order(dims):
    n = len(dims) - 1
    m = [[0] * (n + 1) for _ in range(n + 1)]; s = [[0] * (n + 1) for _ in range(n + 1)]
    for l in range(2, n + 1):
        for i in range(1, n - l + 2):
            j = i + l - 1; m[i][j] = INF
            for k in range(i, j):
                q = m[i][k] + m[k + 1][j] + dims[i - 1] * dims[k] * dims[j]
                if q < m[i][j]: m[i][j], s[i][j] = q, k
    return m, s
def parens(s, i, j, names=None):
    if i == j: return names[i - 1] if names else f"A{i}"
    return "(" + parens(s, i, s[i][j], names) + parens(s, s[i][j] + 1, j, names) + ")"
m, s = matrix_chain_order([30, 35, 15, 5, 10, 20, 25])
print("\nm[1,6] =", m[1][6], " m[2,5] =", m[2][5], " 最优括号化:", parens(s, 1, 6))

# ---------- LCS:原书图 15.8 ----------
def lcs(X, Y):
    m_, n_ = len(X), len(Y); c = np.zeros((m_ + 1, n_ + 1), dtype=int)
    for i in range(1, m_ + 1):
        for j in range(1, n_ + 1):
            c[i, j] = c[i-1, j-1] + 1 if X[i-1] == Y[j-1] else max(c[i-1, j], c[i, j-1])
    out, i, j = [], m_, n_                              # 不用 b 表,直接由 c 回溯(习题 15.4-2)
    while i and j:
        if X[i-1] == Y[j-1]: out.append(X[i-1]); i -= 1; j -= 1
        elif c[i-1, j] >= c[i, j-1]: i -= 1
        else: j -= 1
    return int(c[m_, n_]), "".join(reversed(out))
print("LCS(ABCBDAB, BDCABA) =", lcs("ABCBDAB", "BDCABA"))

# ---------- 最优二叉搜索树:原书图 15.9/15.10 ----------
def optimal_bst(pk, qk):
    n = len(pk) - 1
    e = [[0.0] * (n + 1) for _ in range(n + 2)]; w = [[0.0] * (n + 1) for _ in range(n + 2)]
    root = [[0] * (n + 1) for _ in range(n + 1)]
    for i in range(1, n + 2): e[i][i-1] = w[i][i-1] = qk[i-1]
    for l in range(1, n + 1):
        for i in range(1, n - l + 2):
            j = i + l - 1; e[i][j] = INF; w[i][j] = w[i][j-1] + pk[j] + qk[j]
            for r_ in range(i, j + 1):
                t = e[i][r_-1] + e[r_+1][j] + w[i][j]
                if t < e[i][j]: e[i][j], root[i][j] = t, r_
    return e, root
e, root = optimal_bst([0, .15, .10, .05, .10, .20], [.05, .10, .05, .05, .05, .10])
print(f"最优 BST 期望搜索代价 e[1,5] = {e[1][5]:.2f}; root[1,5]=k{root[1][5]}, root[2,5]=k{root[2][5]}, root[3,5]=k{root[3][5]}")

运行输出:

r = [0, 1, 5, 8, 10, 13, 17, 18, 22, 25, 30]
s = [0, 1, 2, 3, 2, 2, 6, 1, 2, 3, 10]
长度 7 的最优切法: [1, 6]
n=16: 朴素递归  0.007s | 备忘  0.014ms | 自底向上  0.005ms | 三者相等 True
n=20: 朴素递归  0.115s | 备忘  0.023ms | 自底向上  0.010ms | 三者相等 True
n=22: 朴素递归  0.459s | 备忘  0.027ms | 自底向上  0.011ms | 三者相等 True

m[1,6] = 15125  m[2,5] = 7125  最优括号化: ((A1(A2A3))((A4A5)A6))
LCS(ABCBDAB, BDCABA) = (4, 'BCBA')
最优 BST 期望搜索代价 e[1,5] = 2.75; root[1,5]=k2, root[2,5]=k4, root[3,5]=k5

所有数字与原书一致。朴素递归从 \(n=20\) 到 \(n=22\) 时间放大约 4 倍,正是 \(2^n\) 的增长;两种 DP 实现都在微秒级,自底向上比备忘快约一倍(没有递归和散列开销)。

15.8.2 风险模型里的矩阵链:先算什么

因子风险模型把协方差写成 \(\Sigma=BFB^\top+D\),其中 \(B\) 是 \(N\times K\) 的因子暴露、\(F\) 是 \(K\times K\) 的因子协方差。计算 \(M\) 个组合(权重矩阵 \(W\),\(N\times M\))的因子风险 \(W^\top BFB^\top W\) 时,若先把 \(N\times N\) 的 \(BFB^\top\) 算出来,代价和内存都是 \(O(N^2)\) 量级;矩阵链 DP 会告诉你应该先把 \(W\) 投影到因子空间:

import numpy as np, time
INF = float("inf")
def matrix_chain_order(dims):
    n = len(dims) - 1
    m = [[0] * (n + 1) for _ in range(n + 1)]; s = [[0] * (n + 1) for _ in range(n + 1)]
    for l in range(2, n + 1):
        for i in range(1, n - l + 2):
            j = i + l - 1; m[i][j] = INF
            for k in range(i, j):
                q = m[i][k] + m[k + 1][j] + dims[i - 1] * dims[k] * dims[j]
                if q < m[i][j]: m[i][j], s[i][j] = q, k
    return m, s
def parens(s, i, j, names):
    if i == j: return names[i - 1]
    return "(" + parens(s, i, s[i][j], names) + " " + parens(s, s[i][j] + 1, j, names) + ")"
def chain_cost(order, dims):                     # 某个给定括号化的标量乘法次数
    if isinstance(order, int): return 0, (dims[order], dims[order + 1])
    (c1, (a, b)), (c2, (_, d)) = chain_cost(order[0], dims), chain_cost(order[1], dims)
    return c1 + c2 + a * b * d, (a, d)

# 因子风险模型:组合方差 w' B F B' w,N 只股票、K 个因子、M 个组合
N, K, M = 4000, 40, 50
rng = np.random.default_rng(0)
W = rng.normal(size=(N, M)) / N; B = rng.normal(size=(N, K))
L = rng.normal(size=(K, K)); F = L @ L.T / K
names = ["W'", "B", "F", "B'", "W"]; dims = [M, N, K, K, N, M]
m, s = matrix_chain_order(dims)
print("DP 给出的最优顺序:", parens(s, 1, 5, names), f" 标量乘法 {m[1][5]:,}")
naive = (0, ((1, (2, 3)), 4))                    # 先算 N×N 的 B F B',再左右各乘
print("先算 BFB' 的顺序:  ", "(W' ((B (F B')) W))", f" 标量乘法 {chain_cost(naive, dims)[0]:,}")

t0 = time.perf_counter(); V1 = W.T @ ((B @ (F @ B.T)) @ W); t1 = time.perf_counter()
X = B.T @ W; V2 = X.T @ (F @ X); t2 = time.perf_counter()     # 即 (W'B) F (B'W)
print(f"实测:先算 BFB' {t1 - t0:.3f}s;按 DP 顺序 {1e3 * (t2 - t1):.2f}ms;结果一致: {np.allclose(V1, V2)}")

运行输出:

DP 给出的最优顺序: ((W' B) (F (B' W)))  标量乘法 16,180,000
先算 BFB' 的顺序:   (W' ((B (F B')) W))  标量乘法 1,456,400,000
实测:先算 BFB' 0.016s;按 DP 顺序 0.08ms;结果一致: True

推导拆解:手算两种顺序的代价,就能看出差距从哪来(\(N=4000\),\(K=40\),\(M=50\))。 先算 \(BFB^\top\):\(F B^\top\) 是 \(K\times K\) 乘 \(K\times N\),\(40\cdot40\cdot4000=6.4\) 百万;\(B(FB^\top)\) 是 \(N\times K\) 乘 \(K\times N\),\(4000\cdot40\cdot4000=640\) 百万,得到一个 \(4000\times4000\) 的大矩阵;再右乘 \(W\):\(4000\cdot4000\cdot50=800\) 百万;左乘 \(W^\top\):\(50\cdot4000\cdot50=10\) 百万。合计约 1456 百万。 先投影:\(W^\top B\) 是 \(50\cdot4000\cdot40=8\) 百万,得 \(50\times40\);\(B^\top W\) 同样 8 百万,得 \(40\times50\);中间 \(F(B^\top W)\) 只要 \(40\cdot40\cdot50=0.08\) 百万,最后 \((W^\top B)(\cdot)\) 再 0.1 百万。合计约 16.2 百万。 差距的来源一目了然:只要中间结果出现 \(N\times N\),代价就带上 \(N^2\);先把 \(N\) 维的股票空间压到 \(K\) 维的因子空间,所有中间结果都很小。这与你计算组合因子暴露 \(B^\top w\)、再用因子协方差算风险的标准做法完全一致。

标量乘法次数相差约 90 倍,实测时间相差两个数量级(BLAS 多线程让大矩阵乘法有更高的吞吐,所以实测比值与理论比值不完全相同)。更重要的是内存:先算 \(BFB^\top\) 要分配一个 \(4000\times4000\) 的稠密矩阵(128 MB),全市场 5000 只股票、每天几百个组合时,这一步就成了瓶颈。组合优化器里计算风险、风险梯度 \(2\Sigma w\)、边际风险贡献时,都应遵循"先投影到因子空间"的顺序。组合优化的建模与求解见第 04 册,因子风险模型的实战搭建见第 11 册。


本章小结

动态规划把问题分解为相互重叠的子问题,每个子问题只求解一次。适用的两个要素是最优子结构(用剪切—粘贴论证,且子问题必须独立、不争用资源)和重叠子问题(不同子问题个数为多项式)。实现上有带备忘的自顶向下和自底向上两种,渐近复杂度相同;前者只解需要的子问题,后者常数更小、便于压缩空间。运行时间大致等于子问题数乘以每个子问题的选择数。钢条切割、矩阵链乘法、LCS、最优 BST 分别展示了"一维前缀""区间""二维前缀""区间 + 权重"四种典型的子问题空间。

问题 子问题 递归式 复杂度
钢条切割 长度 \(j\) \(r_j=\max_{1\le i\le j}(p_i+r_{j-i})\) \(\Theta(n^2)\)
矩阵链乘法 子链 \(i..j\) \(m[i,j]=\min_k\{m[i,k]+m[k+1,j]+p_{i-1}p_kp_j\}\) \(\Theta(n^3)\) 时间,\(\Theta(n^2)\) 空间
LCS 前缀对 \((i,j)\) 相等取 \(c[i-1,j-1]+1\),否则 \(\max(c[i-1,j],c[i,j-1])\) \(\Theta(mn)\);只求长度 \(O(\min(m,n))\) 空间
最优 BST 关键字区间 \(i..j\) \(e[i,j]=\min_r\{e[i,r-1]+e[r+1,j]+w(i,j)\}\) \(\Theta(n^3)\);Knuth 优化 \(\Theta(n^2)\)
概念 要点
四步法 刻画结构 → 递归定义 → 计算最优值 → 构造最优解
最优子结构 剪切—粘贴;子问题须独立(最长简单路径是反例)
重叠子问题 朴素递归指数时间(钢条 \(2^n\),矩阵链 \(\ge2^{n-1}\))
子问题图 自底向上 = 逆拓扑序;备忘 = 深度优先搜索;时间 ∝ 顶点 + 边

练习

基础

  1. 由 \(T(0)=1\)、\(T(n)=1+\sum_{j=0}^{n-1}T(j)\) 证明 \(T(n)=2^n\)。(原书 15.1-1。)
  2. 举一个反例,说明钢条切割中"每次切单位价格最高的一段"的贪心策略不是最优的。(原书 15.1-2。)
  3. 每次切割有固定成本 \(c\),修改钢条切割的 DP。(原书 15.1-3。)
  4. 对维数序列 \(\langle5,10,3,12,5,50,6\rangle\) 求矩阵链的最优括号化。(原书 15.2-1。)
  5. 求 \(\langle1,0,0,1,0,1,0,1\rangle\) 与 \(\langle0,1,0,1,1,0,1,1,0\rangle\) 的 LCS。(原书 15.4-1。)

进阶

  1. 证明矩阵链长度为 \(n\) 的子问题图有 \(\Theta(n^2)\) 个顶点、\(\Theta(n^3)\) 条边。(原书 15.2-4。)
  2. 举反例说明"先按使 \(p_{i-1}p_kp_j\) 最小的 \(k\) 分割"的贪心不是最优的。(原书 15.3-4。)
  3. 只用 \(\min(m,n)\) 个表项加 \(O(1)\) 额外空间求 LCS 长度。(原书 15.4-4。)
  4. 设计 \(O(n\lg n)\) 的最长递增子序列算法。在一段价格序列上,LIS 的长度能刻画什么?(原书 15.4-6。)
  5. 证明 15.3-5:若每种长度的钢条最多只能卖 \(l_i\) 根,最优子结构不再成立;再想一想,把"剩余配额向量"放进状态能否恢复 DP,代价是什么?
  6. 思考题 15-4(整齐打印):把它改写为"把一段收益率序列切成若干段、使各段内方差之和加分段惩罚最小"的最优分段问题,写出递归式和复杂度。

原书推荐习题:15.1-3、15.1-5、15.2-1、15.2-5、15.3-5、15.3-6、15.4-4、15.4-6、15.5-4,思考题 15-4、15-5、15-7、15-10、15-11。


原书对照

本章小节 原书章节 PDF 页码
15.1 引言、Part IV 导言 Part IV Introduction;Ch.15 开头 p.377–380
15.2 钢条切割 15.1 Rod cutting p.380–391
15.3 矩阵链乘法 15.2 Matrix-chain multiplication p.391–399
15.4 动态规划原理 15.3 Elements of dynamic programming p.399–411
15.5 最长公共子序列 15.4 Longest common subsequence p.411–418
15.6 最优二叉搜索树 15.5 Optimal binary search trees p.418–425
15.7 思考题速览 Problems 15-1~15-12 p.425–433