第 04a 章 分治算法：最大子数组与 Strassen 矩阵乘法

原书第 4 章篇幅较大，本册拆成两章。本章（04a）讲原书 4.1–4.2 节的两个分治算法：一是从"事后最优的一次买卖"引出的最大子数组问题，它与回测报告里的最大回撤是同一个问题的两面；二是 Strassen 矩阵乘法，它说明分治的威力来自"少做一次递归调用"，也说明渐近更快的算法在实践中未必更快。下一章（04b）讲如何求解这些算法产生的递归式，包括代入法、递归树和主定理。

学习目标

1. 把"最大化一次买卖利润"转化为价格变化序列上的最大子数组问题，并说明为什么不能简单地"最低买、最高卖"。
2. 写出最大子数组的分治算法，理解"跨越中点"子问题为何归入合并步骤，得到 \(T(n)=2T(n/2)+\Theta(n)=\Theta(n\lg n)\)。
3. 掌握线性时间的 Kadane 扫描，并把它改写为最大回撤、事后最佳持有区间的计算。
4. 理解朴素分治矩阵乘法为何仍是 \(\Theta(n^3)\)，Strassen 如何用 7 次乘法换来 \(\Theta(n^{\lg7})\)。
5. 说清 Strassen 在实践中常不被采用的四个原因，并能判断量化系统中矩阵运算的真实成本。

读前导读

这一章在解决什么问题

这一章用两个例子把"分治"练熟，并引出一个关键问题：拆成几个子问题、每个多大、合并花多少，怎样决定总耗时。

第一个例子你其实天天见：事后看，一段行情里哪一次买卖最赚钱？ 反过来问，哪一段跌得最狠？ 后者就是回测报告里的最大回撤。最笨的做法是把所有"买入日、卖出日"组合都试一遍，\(n\) 天有约 \(n^2/2\) 种组合，数据翻倍耗时翻 4 倍。分治法把它降到 \(n\lg n\)（翻倍时略多于 2 倍），而一个聪明的线性扫描（Kadane 算法）能降到 \(n\)（翻倍时正好 2 倍）。你在做绩效分析时用"维护历史最高净值、逐日计算回撤"的方法，其实就是这个线性扫描。

第二个例子是矩阵乘法。算协方差矩阵 \(X^\top X\)、因子暴露、组合优化，底层都是矩阵乘法。按定义算，\(n\times n\) 矩阵相乘要 \(n^3\) 次乘法，矩阵边长翻倍耗时翻 8 倍。Strassen 发现一个巧妙的拆法：把矩阵切成四块后，只需 7 次而不是 8 次"半尺寸"乘法，于是翻倍时耗时翻 7 倍，对应 \(n^{\lg 7}\approx n^{2.807}\)。这一章最值得带走的道理是：在递归里，子问题个数比"每次合并多花一点"重要得多；但理论上更快，实践中未必更快。

需要先想起来的数学

• 递归式的读法。 \(T(n)=aT(n/b)+f(n)\) 读作："规模 \(n\) 的耗时 = \(a\) 个规模 \(n/b\) 子问题的耗时 + 拆分与合并的耗时 \(f(n)\)"。例如 \(T(n)=7T(n/2)+\Theta(n^2)\)：把边长 \(n\) 的矩阵乘法变成 7 个边长 \(n/2\) 的乘法，外加一些 \(n^2\) 量级的矩阵加减。怎么解这类式子放在第 04b 章，这一章先用结论。
• 矩阵乘法的定义。 \(C=AB\) 的第 \(i\) 行第 \(j\) 列元素 \(c_{ij}=\sum_k a_{ik}b_{kj}\)，即"\(A\) 的第 \(i\) 行"与"\(B\) 的第 \(j\) 列"对应相乘再相加。例：\(\begin{pmatrix}1&3\\7&5\end{pmatrix}\begin{pmatrix}6&8\\4&2\end{pmatrix}\) 的左上角 \(=1\times6+3\times4=18\)。算一个元素要 \(n\) 次乘法，共 \(n^2\) 个元素，所以是 \(n^3\)。详见 第 00 册第 06 章 线性代数速成。
• 分块矩阵。 把大矩阵切成几块小矩阵后，可以把每一块当成一个"数"来做乘法，规则和 \(2\times2\) 矩阵一样，只是乘法换成矩阵乘法、并且不能交换顺序（\(A_{11}B_{12}\) 不等于 \(B_{12}A_{11}\)）。详见 第 00 册第 06 章 线性代数速成。
• 对数收益可加。 \(\ln\frac{P_j}{P_i}=\ln P_j-\ln P_i=\sum_{t=i+1}^{j}r_t\)，多日对数收益等于每日对数收益之和；简单收益要连乘，不能直接相加。详见 第 00 册第 04 章 级数与收敛 中 e 与指数对数一节。
• \(\lg 7\) 是什么数。 \(\lg 7=\log_27\approx2.807\)，就是"2 的多少次方等于 7"。\(n^{\lg 7}\) 的意思是：\(n\) 翻倍，结果乘 \(2^{\lg 7}=7\)。

怎么读这一章

4.1 节是核心：先看懂"价格问题变成价格变化之和问题"这一步，再看分治的"三种位置"和跨中点子问题，最后看 Kadane 扫描，后者实用价值最大。4.2 节 Strassen 方法的 10 个 \(S\)、7 个 \(P\) 不必记，第一次读只要明白"8 次乘法变成 7 次，为什么指数会下降"，以及"实践中为什么常不用它"。量化实战一（最大回撤）建议细读，量化实战二的结论（"稠密矩阵交给 BLAS，真正要关心的是 \(N^3\)"）对日常工作最有用。相关练习列表第一次可以跳过。

4.0 分治回顾与递归式

分治法在每层递归有三步：分解、解决、合并（第 01 章 1.5 节）。子问题还大、需要递归求解时称递归情况（recursive case）；子问题小到直接求解时，递归"触底"，进入基本情况（base case）。有时除了与原问题同类的更小子问题，还得解一些和原问题不完全相同的子问题，这部分工作算作合并步骤。本章的最大子数组算法就是这样。

**递归式（recurrence）**是用较小输入上的函数值描述函数自身的方程或不等式。归并排序是 \(T(n)=2T(n/2)+\Theta(n)\)。递归式形式多样：按 2/3 与 1/3 不等分割时是 \(T(n)=T(2n/3)+T(n/3)+\Theta(n)\)；子问题也不必按比例缩小，递归线性查找是 \(T(n)=T(n-1)+\Theta(1)\)。求解方法放在第 04b 章。

4.1 最大子数组问题

问题：事后最优的一次买卖

原书的情境是：你可以投资 Volatile Chemical Corporation 的股票，只能买入一股、一次，之后在某天卖出，买卖都在收盘后进行。作为补偿，你能预知未来价格。目标是利润最大。17 天（第 0–16 天）的价格如下（原书图 4.1）：

直觉做法"最低点买、最高点卖"行不通：最低价 63 出现在第 7 天，晚于最高价 113（第 1 天）。退一步，"要么在全局最低点买、要么在全局最高点卖"也不对。原书图 4.2 的反例：价格 \(10,11,7,10,6\)（第 0–4 天），最大利润 3 来自第 2 天以 7 买入、第 3 天以 10 卖出，而 7 不是全局最低（6），10 也不是全局最高（11）。

暴力解：枚举所有"买入早于卖出"的日期对，共 \(\binom n2=\Theta(n^2)\) 对。

变换为最大子数组

换个角度：不看每天的价格，看每天的价格变化。第 \(i\) 天的变化是第 \(i\) 天与第 \(i-1\) 天收盘价之差，得到数组

第 \(i-1\) 天收盘买入、第 \(j\) 天收盘卖出的利润恰好是 \(A[i]+A[i+1]+\cdots+A[j]\)。于是问题变成：找和最大的非空连续子数组，称为最大子数组（maximum subarray）。本例的最大子数组是 \(A[8..11]\)，和为 \(18+20-7+12=43\)：第 7 天收盘以 63 买入，第 11 天收盘以 106 卖出。

推导拆解：为什么"利润 = 一段变化之和"？因为中间的价格两两抵消（这叫裂项相消）。以第 7 天买、第 11 天卖为例： \(A[8]+A[9]+A[10]+A[11]=(P_8-P_7)+(P_9-P_8)+(P_{10}-P_9)+(P_{11}-P_{10})=P_{11}-P_7=106-63=43\)。 金融上这就是"持有期损益 = 每日盯市损益之和"，和期货逐日结算的道理一样：每天的变动单独记账，加总起来就是总损益。好处是：原来要同时决定两个日期，现在变成在一个数组里找"和最大的一段"，可以用统一的算法处理。

变换本身不降低复杂度：\(n\) 天仍有 \(\Theta(n^2)\) 个子数组。利用已算出的部分和，每个子数组的和可以 \(O(1)\) 得到，暴力法是 \(\Theta(n^2)\)（练习 4.1-2）。最大子数组可能不唯一，所以说"一个"最大子数组。只有数组含负数时问题才有意义，全为非负时整个数组就是答案。

分治解法

求 \(A[low..high]\) 的最大子数组。取中点 \(mid\)，任何连续子数组 \(A[i..j]\) 必定恰好落在三处之一：

1. 完全在左半 \(A[low..mid]\)：\(low\le i\le j\le mid\)；
2. 完全在右半 \(A[mid+1..high]\)：\(mid<i\le j\le high\)；
3. 跨越中点：\(low\le i\le mid<j\le high\)。

前两种是规模减半的同类问题，递归求解。第三种不是原问题的更小实例，因为多了"必须跨过中点"的限制，所以把它归入合并步骤。好在它可以在线性时间内解决：任何跨中点的子数组都由 \(A[i..mid]\) 和 \(A[mid+1..j]\) 拼成，两部分可以独立地各自取最大，再相加。

FIND-MAX-CROSSING-SUBARRAY(A, low, mid, high)
 1  left-sum = -∞
 2  sum = 0
 3  for i = mid downto low
 4      sum = sum + A[i]
 5      if sum > left-sum
 6          left-sum = sum
 7          max-left = i
 8  right-sum = -∞
 9  sum = 0
10  for j = mid + 1 to high
11      sum = sum + A[j]
12      if sum > right-sum
13          right-sum = sum
14          max-right = j
15  return (max-left, max-right, left-sum + right-sum)

左半循环从 \(mid\) 往下走，考虑的每个子数组都形如 \(A[i..mid]\)；右半从 \(mid+1\) 往上走，都形如 \(A[mid+1..j]\)。两个循环共迭代 \((mid-low+1)+(high-mid)=n\) 次，每次常数时间，所以是 \(\Theta(n)\)。

推导拆解：用本节的 16 个价格变化走一遍顶层调用，\(low=1\)、\(high=16\)、\(mid=8\)。 左半从第 8 位往回累加：\(18,\ -5,\ -21,\ -24,\ -4,\ -29,\ -32,\ -19\)，最大的是 18，对应 \(i=8\)。意思是"跨中点的那一段，左边只带上 \(A[8]\) 最划算"。 右半从第 9 位往后累加：\(20,\ 13,\ 25,\ 20,\ -2,\ 13,\ 9,\ 16\)，最大的是 25，对应 \(j=11\)。 跨中点的最佳是 \(18+25=43\)，即 \(A[8..11]\)。递归求得左半 \(A[1..8]\) 的最佳是单独的 \(A[4]=20\)，右半 \(A[9..16]\) 的最佳是 \(A[9..11]=25\)。三者取最大，答案是跨中点的 43。 为什么左右可以各自取最大再相加？因为跨中点的一段一定是"以 \(mid\) 结尾的左段"加"以 \(mid+1\) 开头的右段"，两段的选择互不影响，总和最大就等于各自最大之和。就像一笔跨年度的持仓收益 = 去年那部分 + 今年那部分，两部分可以分开优化。 金融直觉：在价格图上，\(mid\) 那天的收盘价是个固定的"中转站"。左段问"在 \(mid\) 之前哪天买入，到 \(mid\) 时赚得最多"，右段问"从 \(mid\) 起持有到哪天卖出赚得最多"。

FIND-MAXIMUM-SUBARRAY(A, low, high)
 1  if high == low
 2      return (low, high, A[low])              // 基本情况：只有一个元素
 3  else mid = ⌊(low + high)/2⌋
 4      (left-low, left-high, left-sum) = FIND-MAXIMUM-SUBARRAY(A, low, mid)
 5      (right-low, right-high, right-sum) = FIND-MAXIMUM-SUBARRAY(A, mid + 1, high)
 6      (cross-low, cross-high, cross-sum) = FIND-MAX-CROSSING-SUBARRAY(A, low, mid, high)
 7      if left-sum >= right-sum and left-sum >= cross-sum
 8          return (left-low, left-high, left-sum)
 9      elseif right-sum >= left-sum and right-sum >= cross-sum
10          return (right-low, right-high, right-sum)
11      else return (cross-low, cross-high, cross-sum)

分析。 设 \(n\) 为 2 的幂。基本情况 \(T(1)=\Theta(1)\)。递归情况：第 1、3 行常数；两个规模 \(n/2\) 的子问题共 \(2T(n/2)\)；跨中点 \(\Theta(n)\)；第 7–11 行常数。于是

与归并排序的递归式完全相同，解为 \(\Theta(n\lg n)\)，递归栈空间 \(\Theta(\lg n)\)。分治法渐近快于暴力法，但它不是最优的。

线性时间：Kadane 扫描（练习 4.1-5）

从左往右扫描，同时记住两样东西：目前为止的最大子数组，以及以当前位置结尾的最大子数组。关键观察是：以 \(j+1\) 结尾的最大子数组，要么只含 \(A[j+1]\)，要么是"以 \(j\) 结尾的最大子数组"再接上 \(A[j+1]\)。若以 \(j\) 结尾的最佳和不大于 0，接上它只会更差，就从 \(j+1\) 重新开始。每步常数时间，总计 \(\Theta(n)\)，额外空间 \(O(1)\)。这是一个动态规划的雏形（原书第 15 章）。

金融直觉：\(\text{cur}_j\) 可以理解为"如果今天必须卖出，最好的入场点能带来的累计盈利"。每天只做一个决定：昨天的累计盈利如果已经是负的（或零），说明之前任何时点入场到昨天都是亏的，那不如忘掉过去、按昨天收盘价重新入场，所以今天的累计只算今天的涨跌 \(A[j+1]\)；如果昨天的累计是正的，就把今天的涨跌接上去。\(\text{best}\) 记录扫描过程中见过的最好成绩。 用前几天走一遍：\(A=\langle13,-3,-25,20,\dots\rangle\)。第 1 天 cur \(=13\)；第 2 天 \(13-3=10\)；第 3 天 \(10-25=-15\)；第 4 天，昨天累计为负，重新开始，cur \(=20\)。best 依次是 \(13,13,13,20\)。 它比分治快的原因是：分治在每层都重新扫一遍跨中点的情况，而 Kadane 把"以每一天结尾的最佳结果"记下来，后一天直接复用前一天的结果，没有重复劳动。

其他几道练习的结论：全为负数时，最大子数组就是值最大的那个单元素（4.1-1）；若允许空子数组（和为 0），只需在结果为负时改为返回空子数组（4.1-4）；练习 4.1-3 让读者实测暴力法与分治法的交叉点 \(n_0\)，再把分治的基本情况改为 \(n<n_0\) 时调用暴力法，这和第 01 章的"粗化叶子"是同一个思想。

4.2 Strassen 矩阵乘法

直接法与朴素分治

\(n\times n\) 矩阵 \(A=(a_{ij})\)、\(B=(b_{ij})\) 的乘积 \(C=AB\) 的元素是

三重循环各 \(n\) 次，直接法 \(\Theta(n^3)\)。直觉上似乎任何矩阵乘法都要 \(\Omega(n^3)\)，但这是错的。

先看一个简单的分治。设 \(n\) 是 2 的幂，把 \(A\)、\(B\)、\(C\) 各分成四个 \(n/2\times n/2\) 块：

划分可以用下标计算在 \(\Theta(1)\) 内完成（用原矩阵的行列下标范围表示子矩阵，就像 NumPy 的切片视图；即使复制成新矩阵也只是 \(\Theta(n^2)\)，不改变结论）。之后做 8 次规模 \(n/2\) 的递归乘法和 4 次 \(n/2\times n/2\) 矩阵加法（每次 \(\Theta(n^2/4)\)）：

解为 \(\Theta(n^3)\)，并不比直接法快。

白话解释：分块乘法的公式和 \(2\times2\) 数字矩阵的乘法一模一样：\(C_{11}=A_{11}B_{11}+A_{12}B_{21}\) 就是"\(A\) 的第一行块"乘"\(B\) 的第一列块"。唯一的区别是每个字母现在代表一个 \(n/2\times n/2\) 的小矩阵，乘号是矩阵乘法，所以顺序不能颠倒。四个 \(C\) 块各需要 2 次小乘法，共 8 次。可以用 \(n=2\) 自检：每块是 \(1\times1\)，公式退化成普通的 \(c_{11}=a_{11}b_{11}+a_{12}b_{21}\)。

一个容易犯的错误。 \(\Theta\) 记号会吸收常数，划分的 \(\Theta(2)\)、加法的 \(\Theta(4\cdot n^2/4)\) 都写成 \(\Theta(1)\)、\(\Theta(n^2)\)。但 \(T(n/2)\) 前面的系数 8 不能吸收：它决定递归树每个结点有几个孩子，从而决定每层有多少项。去掉它，递归树就从"茂密"变成一条链，结论完全不同。

推导拆解：用"边长翻倍，耗时翻几倍"来看系数 8 和 7 的作用。 边长减半，一次小乘法的工作量是原来的多少？直接法的工作量约 \(n^3\)，边长减半后是 \((n/2)^3=n^3/8\)。所以 8 个半尺寸乘法加起来 \(8\times n^3/8=n^3\)，一点没省，这就是朴素分治仍是 \(\Theta(n^3)\) 的原因：拆完以后，总工作量原封不动。 下面要讲的 Strassen 只做 7 个：每往下拆一层，乘法工作量就乘 \(7/8\)。拆 \(\lg n\) 层一直拆到 \(1\times1\)，总共做 \(7^{\lg n}\) 次标量乘法。用 3.3 节的恒等式 \(a^{\log_bc}=c^{\log_ba}\)，\(7^{\lg n}=n^{\lg 7}\approx n^{2.807}\)。翻倍时耗时乘 7 而不是 8。 它多出来的 18 次矩阵加减（构造 10 个 \(S\)、组合 \(C\) 的 8 次）呢？每次加减是 \(\Theta(n^2)\)，边长翻倍只翻 4 倍，比 7 倍慢，所以在规模大时被乘法部分压住，不影响指数。这就是"系数 \(a\) 不能丢、\(f(n)\) 里的常数可以丢"的直观原因：\(a\) 是每一层都会复利一次的倍数，\(f(n)\) 的常数只是一次性的比例。

Strassen 方法

关键想法是让递归树稍微不那么茂密：只做 7 次而不是 8 次 \(n/2\times n/2\) 的递归乘法，代价是多做常数次矩阵加减法。四步：

1. 把 \(A\)、\(B\)、\(C\) 划分为 \(n/2\times n/2\) 子矩阵，\(\Theta(1)\)。
2. 构造 10 个矩阵，每个是两个子矩阵的和或差，共 \(\Theta(n^2)\)：

3. 递归计算 7 个乘积：

（每行等号右边的展开只是为了验证，真正要算的只有中间那一次乘法。）

4. 用 \(P_i\) 的加减得到 \(C\) 的四块（共 8 次加减，\(\Theta(n^2)\)）：

验证 \(C_{11}\)：把 \(P_5+P_4-P_2+P_6\) 展开，\(A_{11}B_{22}\)、\(A_{22}B_{11}\)、\(A_{22}B_{22}\)、\(A_{22}B_{21}\)、\(A_{12}B_{22}\) 各出现一正一负，相互抵消，剩下 \(A_{11}B_{11}+A_{12}B_{21}\)。其余三块类似。

白话解释：Strassen 的思路是"用便宜的加法换贵的乘法"。先看一个小版本，也就是练习 4.2-7 的复数乘法：\((a+bi)(c+di)\) 按定义要 4 次乘法 \(ac,bd,ad,bc\)。但虚部 \(ad+bc\) 可以写成 \((a+b)(c+d)-ac-bd\)，而 \(ac\)、\(bd\) 反正要算，于是只需 3 次乘法，代价是多做几次加减。Strassen 对 \(2\times2\) 分块矩阵做了同样的事：精心挑选 7 个"先加减、再相乘"的组合 \(P_1,\dots,P_7\)，使得 4 个 \(C\) 块都能由它们加减得出。 这在数字世界里不值得（乘一个数和加一个数差不多快），但当每个"数"是一个大矩阵时，一次矩阵乘法约 \(n^3\)、一次矩阵加法约 \(n^2\)，用十几次加法换掉一次乘法就很划算，而且这个交换在递归的每一层都会发生。 再用练习 4.2-1 的数字检查 \(C_{12}=P_1+P_2\)：\(P_1=A_{11}(B_{12}-B_{22})=1\times(8-2)=6\)，\(P_2=(A_{11}+A_{12})B_{22}=(1+3)\times2=8\)，和为 14，正是答案矩阵右上角。

递归式变为

由主定理（第 04b 章）得 \(T(n)=\Theta(n^{\lg7})\)，\(\lg7\approx2.807\)。额外空间每层 \(\Theta(n^2)\) 临时矩阵，几何级数求和后总计 \(\Theta(n^2)\)。原书作者说 Strassen 的构造"一点也不显然"，并称这可能是全书最大的轻描淡写。

相关练习

• 练习 4.2-1：用 Strassen 算 \(\begin{pmatrix}1&3\\7&5\end{pmatrix}\begin{pmatrix}6&8\\4&2\end{pmatrix}=\begin{pmatrix}18&14\\62&66\end{pmatrix}\)（下面的代码会验证）。
• 练习 4.2-3：\(n\) 不是 2 的幂时，补零扩充到不小于 \(n\) 的 2 的幂。尺寸至多翻倍，仍是 \(\Theta(n^{\lg7})\)。
• 练习 4.2-4：若能用 \(k\) 次乘法（不假设乘法交换）完成 \(3\times3\) 矩阵乘法，递归后为 \(\Theta(n^{\log_3k})\)。要比 Strassen 快需 \(\log_3k<\lg7\)，即 \(k<3^{\lg 7}\approx21.85\)，最大 \(k=21\)，此时约为 \(n^{2.77}\)。
• 练习 4.2-5：V. Pan 的方法用 132,464 次乘法做 \(68\times68\)，143,640 次做 \(70\times70\)，155,424 次做 \(72\times72\)。指数分别是 \(\log_{68}132464\approx2.795128\)、\(\log_{70}143640\approx2.795123\)、\(\log_{72}155424\approx2.795147\)，\(70\times70\) 最好，三者都优于 Strassen 的 2.807。
• 练习 4.2-6：用 Strassen 作子程序，\(kn\times n\) 乘 \(n\times kn\) 需 \(k^2\) 次 \(n\times n\) 乘法，\(\Theta(k^2n^{\lg7})\)；反过来 \(n\times kn\) 乘 \(kn\times n\) 只需 \(k\) 次乘法再求和，\(\Theta(kn^{\lg7})\)。
• 练习 4.2-7：只用三次实数乘法计算复数乘积 \((a+bi)(c+di)\)：算 \(ac\)、\(bd\)、\((a+b)(c+d)\)，实部 \(ac-bd\)，虚部 \((a+b)(c+d)-ac-bd\)。这和 Karatsuba 大整数乘法是同一个技巧。

Strassen 在实践中的位置

Strassen 1969 年发表时引起轰动，此前很少有人想到存在渐近快于直接法的矩阵乘法。之后上界不断改进，原书写作时最好的是 Coppersmith–Winograd 的 \(O(n^{2.376})\)；已知最好的下界只是显然的 \(\Omega(n^2)\)（至少要填满 \(n^2\) 个元素）。

但实践中 Strassen 常不是首选，原书章末注记列了四个原因：

1. \(\Theta(n^{\lg7})\) 中隐藏的常数因子比直接法大；
2. 稀疏矩阵有专门的更快方法；
3. 数值稳定性不如直接法，浮点运算误差积累更大；
4. 各层递归的临时子矩阵占用空间。

后两点在 1990 年前后有所缓解：Higham 指出稳定性的差别被夸大了，对某些应用不可接受、对另一些可以接受；Bailey、Lee、Simon 讨论了降低内存的技术。实际的稠密矩阵乘法实现会在规模高于某个交叉点时用 Strassen，低于时切回直接法。交叉点高度依赖系统：只数运算次数、忽略 cache 与流水线的分析给出 \(n=8\)（Higham）或 \(n=12\)（Huss-Lederman 等）；D'Alberto 与 Nicolau 的自适应方案在安装时实测，不同系统上的交叉点从 \(n=400\) 到 \(2150\) 不等，有的系统甚至找不到交叉点。

量化实战一：最大回撤与事后最佳持有区间

最大回撤就是"最小子数组"

**最大回撤（maximum drawdown）**是净值从某个历史高点到其后最低点的最大跌幅，是回测报告的标准风险指标。设 \(p_t=\ln P_t\) 为对数净值，\(r_t=p_t-p_{t-1}\) 为对数收益。从第 \(i\) 天到第 \(j\) 天的对数涨跌是 \(r_{i+1}+\cdots+r_j\)，所以：

• 事后最佳的一次持有区间 = 对数收益序列的最大子数组；
• 最大回撤（对数尺度）= 对数收益序列的最小子数组和 \(m\)，换算成百分比为 \(1-e^{m}\)。

这里必须用对数收益，因为对数收益可以相加，简单收益不行。最小子数组等于对 \(-r\) 求最大子数组再取负。实务中更常用的写法是维护"历史最高点"，每天计算当前回撤，这与 Kadane 扫描本质相同，都是 \(\Theta(n)\)。

要强调的是，最大回撤和事后最佳区间都是**事后（look-ahead）**概念，只能用来评估一段已经发生的历史，不能作为交易信号。原书假设"能预知未来价格"，正是为了把这一点说清楚。

推导拆解：为什么回撤百分比是 \(1-e^{m}\)？设峰值日净值 \(P_i\)、谷底日 \(P_j\)，对数尺度上的跌幅是 \(m=\ln P_j-\ln P_i=r_{i+1}+\cdots+r_j\)（一段对数收益之和，为负）。回撤百分比 \(=1-\frac{P_j}{P_i}=1-e^{\ln P_j-\ln P_i}=1-e^{m}\)。例：\(m=-0.693\) 时 \(e^m=0.5\)，回撤 50%。最大回撤要找"和最小的一段"，就是最小子数组。 简单收益为什么不行？净值从 100 跌 50% 到 50、再涨 50% 到 75，简单收益相加是 0，实际却亏了 25%；对数收益 \(\ln0.5+\ln1.5=\ln0.75\)，加起来正好是真实结果。

代码说明：np.diff(price) 计算相邻两天的差，得到价格变化数组；range(mid, lo - 1, -1) 表示从 mid 倒数到 lo（第三个参数 -1 是步长）；max(a, b, c, key=lambda t: t[0]) 在三个结果里按"每个结果的第 0 项（和）"比大小，取最大的那个；np.cumsum(r) 是累计求和，把对数收益变成对数净值。max_drawdown 函数用的正是"维护历史高点 peak"的写法。

import time
import numpy as np

def max_sub_brute(A):                      # Θ(n^2)，练习 4.1-2
    best = (-np.inf, 0, 0)
    for i in range(len(A)):
        s = 0.0
        for j in range(i, len(A)):
            s += A[j]
            if s > best[0]:
                best = (s, i, j)
    return best

def max_crossing(A, lo, mid, hi):          # 跨中点的最大子数组，Θ(n)
    ls, s, ml = -np.inf, 0.0, mid
    for i in range(mid, lo - 1, -1):
        s += A[i]
        if s > ls: ls, ml = s, i
    rs, s, mr = -np.inf, 0.0, mid + 1
    for j in range(mid + 1, hi + 1):
        s += A[j]
        if s > rs: rs, mr = s, j
    return (ls + rs, ml, mr)

def max_sub_dc(A, lo, hi):                 # 分治，T(n)=2T(n/2)+Θ(n)=Θ(n lg n)
    if lo == hi:
        return (A[lo], lo, hi)
    mid = (lo + hi) // 2
    return max(max_sub_dc(A, lo, mid), max_sub_dc(A, mid + 1, hi),
               max_crossing(A, lo, mid, hi), key=lambda t: t[0])

def max_sub_kadane(A):                     # 线性扫描，Θ(n)，练习 4.1-5
    best, cur, start = (-np.inf, 0, 0), 0.0, 0
    for j, a in enumerate(A):
        if cur <= 0:
            cur, start = a, j
        else:
            cur += a
        if cur > best[0]:
            best = (cur, start, j)
    return best

price = [100,113,110,85,105,102,86,63,81,101,94,106,101,79,94,90,97]   # 图 4.1
A = np.diff(price).astype(float)           # A[0] 对应原书 A[1]
for f in (max_sub_brute, lambda a: max_sub_dc(a, 0, len(a) - 1), max_sub_kadane):
    s, i, j = f(A)
    print(f"利润 {s:.0f}: 第 {i} 天收盘后买入(价 {price[i]}), 第 {j+1} 天收盘后卖出(价 {price[j+1]})")

def max_drawdown(logp):                    # 最大回撤 = 对数收益的“最小子数组和”
    peak, peak_t, mdd, span = logp[0], 0, 0.0, (0, 0)
    for t in range(1, len(logp)):
        if logp[t] > peak:
            peak, peak_t = logp[t], t
        elif logp[t] - peak < mdd:
            mdd, span = logp[t] - peak, (peak_t, t)
    return 1 - np.exp(mdd), span

rng = np.random.default_rng(7)
r = rng.normal(0.0003, 0.012, 2520)        # 约 10 年日对数收益
logp = np.concatenate([[0.0], np.cumsum(r)])
dd, (p, q) = max_drawdown(logp)
s, i, j = max_sub_kadane(-r)               # 最小子数组和 = -(对 -r 的最大子数组和)
print(f"最大回撤 {dd:.2%}，峰值日 {p}，谷底日 {q}；Kadane(-r) 得 {1-np.exp(-s):.2%}，区间 {i}..{j+1}")
s, i, j = max_sub_kadane(r)
print(f"事后最佳单次持有：第 {i} 日买入、第 {j+1} 日卖出，收益 {np.exp(s)-1:.2%}")

for n in (500, 1000, 2000):
    x = rng.standard_normal(n)
    ts = []
    for f in (max_sub_brute, lambda a: max_sub_dc(a, 0, len(a) - 1), max_sub_kadane):
        t0 = time.perf_counter(); f(x); ts.append(time.perf_counter() - t0)
    print(f"n={n:5d}  暴力 {ts[0]*1e3:7.1f} ms  分治 {ts[1]*1e3:6.2f} ms  Kadane {ts[2]*1e3:5.2f} ms")

输出（计时随机器而变）：

利润 43: 第 7 天收盘后买入(价 63), 第 11 天收盘后卖出(价 106)
利润 43: 第 7 天收盘后买入(价 63), 第 11 天收盘后卖出(价 106)
利润 43: 第 7 天收盘后买入(价 63), 第 11 天收盘后卖出(价 106)
最大回撤 59.60%，峰值日 2，谷底日 854；Kadane(-r) 得 59.60%，区间 2..854
事后最佳单次持有：第 854 日买入、第 2159 日卖出，收益 109.54%
n=  500  暴力     6.1 ms  分治   0.43 ms  Kadane  0.03 ms
n= 1000  暴力    25.1 ms  分治   0.93 ms  Kadane  0.07 ms
n= 2000  暴力   100.5 ms  分治   1.98 ms  Kadane  0.12 ms

三种算法在原书数据上都给出利润 43、第 7 天买入第 11 天卖出。在模拟的十年净值上，"维护历史高点"与"对 \(-r\) 做 Kadane"得到相同的最大回撤和区间。计时显示 \(n\) 翻倍时暴力法约变为 4 倍、分治约 2 倍多一点、Kadane 约 2 倍，分别对应 \(n^2\)、\(n\lg n\)、\(n\)。

还有一点值得注意：这条模拟路径的年化漂移约 7.5%、年化波动约 19%，最大回撤却接近 60%。随机游走路径在十年里出现很深的回撤并不罕见，单看一条回测曲线的回撤会低估风险。

其他量化对应

• 滚动最大回撤：在每个窗口内都重新算是 \(\Theta(nw)\)；用单调队列维护窗口内最大值可以把"窗口内当前回撤"降到 \(\Theta(n)\)。
• 多资产的最佳区间：对 \(m\) 个资产各做一遍 Kadane 是 \(\Theta(mn)\)，向量化后可以一次扫完。

量化实战二：矩阵乘法的真实成本

矩阵乘法是协方差估计、因子暴露计算、组合优化的核心运算。下面的代码验证练习 4.2-1，比较标量乘法次数，并在 \(1024\times1024\) 的收益矩阵上比较 NumPy（调用 BLAS）与纯 NumPy 写的 Strassen。

import time
import numpy as np

def strassen(A, B, leaf=64):
    """Strassen 乘法（n 为 2 的幂）。n <= leaf 时退回普通乘法：粗化递归叶子。"""
    n = A.shape[0]
    if n <= leaf:
        return A @ B
    h = n // 2
    A11, A12, A21, A22 = A[:h, :h], A[:h, h:], A[h:, :h], A[h:, h:]   # 视图，Θ(1) 划分
    B11, B12, B21, B22 = B[:h, :h], B[:h, h:], B[h:, :h], B[h:, h:]
    P1 = strassen(A11, B12 - B22, leaf)
    P2 = strassen(A11 + A12, B22, leaf)
    P3 = strassen(A21 + A22, B11, leaf)
    P4 = strassen(A22, B21 - B11, leaf)
    P5 = strassen(A11 + A22, B11 + B22, leaf)
    P6 = strassen(A12 - A22, B21 + B22, leaf)
    P7 = strassen(A11 - A21, B11 + B12, leaf)
    C = np.empty_like(A)
    C[:h, :h] = P5 + P4 - P2 + P6
    C[:h, h:] = P1 + P2
    C[h:, :h] = P3 + P4
    C[h:, h:] = P5 + P1 - P3 - P7
    return C

# 练习 4.2-1
A = np.array([[1., 3.], [7., 5.]]); B = np.array([[6., 8.], [4., 2.]])
print(strassen(A, B, leaf=1))

# 标量乘法次数：朴素分治 M(n)=8M(n/2)，Strassen M(n)=7M(n/2)，M(1)=1
for k in (4, 7, 10):
    n = 2 ** k
    print(f"n={n:5d}: 普通 {n**3:>13,d} 次, Strassen {7**k:>13,d} 次, 比值 {7**k/n**3:.3f}")

# 收益率矩阵算协方差的场景：X 为 T×N 去均值收益，Σ = X'X/(T-1)
rng = np.random.default_rng(1)
N = 1024
X = rng.standard_normal((N, N))
t0 = time.perf_counter(); S1 = X.T @ X; t1 = time.perf_counter()
S2 = strassen(np.ascontiguousarray(X.T), X, leaf=128); t2 = time.perf_counter()
rel = np.abs(S1 - S2).max() / np.abs(S1).max()
print(f"BLAS {1e3*(t1-t0):.1f} ms, Strassen(numpy) {1e3*(t2-t1):.1f} ms, 最大相对误差 {rel:.1e}")

输出（计时随机器而变）：

[[18. 14.]
 [62. 66.]]
n=   16: 普通         4,096 次, Strassen         2,401 次, 比值 0.586
n=  128: 普通     2,097,152 次, Strassen       823,543 次, 比值 0.393
n= 1024: 普通 1,073,741,824 次, Strassen   282,475,249 次, 比值 0.263
BLAS 5.3 ms, Strassen(numpy) 21.5 ms, 最大相对误差 2.5e-15

白话解释：代码里的 A @ B 是 Python 的矩阵乘法运算符，背后调用 BLAS 这类高度优化的数值库；A[:h, :h] 取左上角 \(h\times h\) 的块，返回的是"视图"（指向原数据的一个窗口，不复制）。乘法次数表里的比值 \(7^k/8^k=(7/8)^k\) 随层数 \(k\) 一路下降：每多拆一层，就省下 1/8 的乘法。

\(n=1024\) 时 Strassen 的标量乘法只有直接法的 26%，但实测反而比 BLAS 慢约 4 倍。原因正是原书列出的几条：大量临时矩阵的分配与加减（内存带宽），以及 BLAS 对 cache 分块、SIMD 向量化、多线程的深度优化，这些常数因子 RAM 模型都看不到。量化开发中的结论是：

• 稠密矩阵运算一律交给 BLAS/LAPACK（NumPy、SciPy 背后就是它们），不要自己写分治乘法；
• 真正要关心的是 \(\Theta(n^3)\) 这个指数本身：几千只股票的协方差求逆、组合优化里的矩阵分解都是 \(N^3\) 量级，股票池扩大 10 倍成本涨 1000 倍。这时的出路是改变问题结构，例如用 \(K\) 个因子的因子模型把协方差写成 \(B\Sigma_fB^\top+D\)，求逆借助 Woodbury 公式降到 \(O(NK^2+K^3)\)，或者利用稀疏性、低秩近似；
• 思考题 4-2 提醒的参数传递问题在 Python 里很常见：NumPy 切片是视图（\(\Theta(1)\)），花式索引和 .copy() 是复制（\(\Theta(n)\)）。上面的 Strassen 用切片视图做划分，符合原书"用下标计算在 \(\Theta(1)\) 内划分"的假设；若在递归里反复复制大数组或按值传递大 DataFrame，原本 \(\Theta(N\lg N)\) 的递归可能退化成 \(\Theta(N^2)\)。

本章小结

分治法的威力取决于两件事：子问题能否独立求解，合并代价有多大。最大子数组问题中，左右两半递归求解，跨越中点的情况归入合并步骤，用两次线性扫描解决，得到与归并排序相同的 \(T(n)=2T(n/2)+\Theta(n)=\Theta(n\lg n)\)。但分治不一定最优，Kadane 扫描只需 \(\Theta(n)\)，并且直接给出了最大回撤和事后最佳持有区间的计算方法（要用对数收益）。矩阵乘法中，朴素分治的 8 次递归调用仍是 \(\Theta(n^3)\)；Strassen 用 7 次乘法加常数次矩阵加减，把指数降到 \(\lg7\approx2.807\)。递归式中子问题个数 \(a\) 不能被 \(\Theta\) 吸收，它决定递归树的分支数。实践中渐近优势要越过交叉点才能兑现，Strassen 因常数、稳定性和内存原因常不被采用，量化系统应依赖 BLAS，并从问题结构上降低 \(N^3\) 的成本。

练习

基础

1. 所有元素都为负数时，FIND-MAXIMUM-SUBARRAY 返回什么？ 答案：值最大（绝对值最小）的那个单元素。
2. 写出最大子数组的 \(\Theta(n^2)\) 暴力算法，要求每个子数组的和 \(O(1)\) 得到。
3. 用 Strassen 算法手算 \(\begin{pmatrix}1&3\\7&5\end{pmatrix}\begin{pmatrix}6&8\\4&2\end{pmatrix}\)，写出 \(S_1,\dots,S_{10}\) 与 \(P_1,\dots,P_7\)。 答案：\(\begin{pmatrix}18&14\\62&66\end{pmatrix}\)。
4. 只用三次实数乘法计算 \((a+bi)(c+di)\)。
5. 修改 Kadane 算法使其允许返回空子数组（和为 0），并说明这对应"整个区间都不交易"。

进阶

6. 实现暴力法与分治法，实测交叉点 \(n_0\)；再把分治的基本情况改为 \(n<n_0\) 时调用暴力法，观察交叉点是否改变（练习 4.1-3）。
7. 给定日对数收益序列，用一次 \(\Theta(n)\) 扫描同时输出：最大回撤及其峰谷日期、回撤后恢复到前高所需天数。 提示：在维护历史高点的同时，记录回撤开始后第一次重新达到该高点的日期。
8. 若能用 \(k\) 次乘法完成 \(3\times3\) 矩阵乘法，求使 \(n\times n\) 乘法为 \(o(n^{\lg7})\) 的最大 \(k\)。 答案：\(k=21\)，运行时间 \(\Theta(n^{\log_321})\approx n^{2.77}\)。
9. 允许多次交易（每次持有不重叠、卖出后才能再买）时，事后最大利润是多少？与本章的单次买卖相比，复杂度如何？ 提示：把所有正的日价格变化相加即可，\(\Theta(n)\)。这说明问题约束的小变化可以让问题变简单。
10. 设 \(X\) 为 \(T\times N\) 收益矩阵。分别估算计算 \(X^\top X\)（\(N\times N\)）与 \(XX^\top\)（\(T\times T\)）的运算量；当 \(N=5000\)、\(T=250\) 时，哪种做法更适合用来求协方差矩阵的特征分解？ 提示：前者 \(\Theta(N^2T)\)、后者 \(\Theta(T^2N)\)；两个矩阵的非零特征值相同，\(T\ll N\) 时对 \(T\times T\) 矩阵做分解再还原特征向量更便宜。

原书推荐习题：4.1-5（线性时间最大子数组，并改写为最大回撤）；4.1-3（交叉点实测）；4.2-7（三次乘法算复数积，与 Karatsuba 相通）；4.2-4、4.2-5（比较递归乘法的指数）；思考题 4-2（参数传递代价）。

原书对照