第 07 章 导数计算

学习目标

读完本章，你应当能够：

1. 推导前向差分与中心差分的截断误差和舍入误差，给出最优步长（\(\sqrt u\) 与 \(u^{1/3}\)）及可达精度，并会在有噪声的函数（如蒙特卡洛定价）上正确使用差分。
2. 用有限差分计算 Jacobian–向量积和 Hessian–向量积，理解它们为什么只需一次额外的函数或梯度求值。
3. 用图着色把稀疏 Jacobian/Hessian 的差分次数降到每行非零元个数的量级，并会利用 Hessian 的对称性进一步节省。
4. 理解自动微分的前向模式与反向模式：计算图、种子向量、伴随变量；说出反向模式"梯度代价至多为函数代价的 4–5 倍、与变量个数无关"这一关键性质及其存储代价。
5. 知道用自动微分计算 Jacobian、Hessian 和 Hessian–向量积的方法与代价，以及自动微分的局限（截断误差的导数、代码分支）。
6. 把本章方法用于期权 Greeks 计算：有限差分步长选择、公共随机数、伴随算法微分（AAD）。

读前导读

这一章在解决什么问题

一句话：怎样让计算机又快又准地算出导数。前几章的优化算法（最速下降、牛顿法、信赖域）都默认"梯度和 Hessian 是现成的"，这一章回答它们从哪里来。

你其实每天都在用导数，只是名字不同。Delta 是期权价格对标的价格的导数，Vega 是对波动率的导数，久期是债券价格对收益率的导数（再除以价格、取负号），DV01 是"收益率动 1bp 价格变多少"。风险系统算 DV01 的标准做法是"把收益率曲线上移 1bp，重新定价一遍，看价格差"——这就是本章的有限差分。它简单，但有两个问题：步长取多大才合适（太大不准，太小被计算机的舍入误差淹没），以及参数一多就很慢（一本有 500 个曲线节点的利率互换账簿，要重新定价 500 次以上）。

本章的第二个主角自动微分（AD）专门解决"参数多"的问题。它把定价程序拆成一串加减乘除、exp、log 等基本运算，对每一步用链式法则，最终得到的是精确导数而不是近似。其中的反向模式有一个惊人的性质：不管有多少个输入参数，算出全部一阶敏感度的代价只是一次定价的几倍。投行做 XVA、大型账簿的风险计算所用的 AAD（伴随算法微分），就是反向模式 AD。

需要先想起来的数学

1. 导数与差分。 导数 \(f'(x)\) 是"自变量动一点点，函数值动多少"的比率极限：\(f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}{h}\)。有限差分就是不取极限、直接用一个小的 \(h\) 算这个比率。例：\(f(x)=x^2\)，\(x=3\)，\(h=0.01\)，差分 \(=(9.0601-9)/0.01=6.01\)，真值 6，误差 0.01 正好约等于 \(h\)。见 第 00 册第 02 章 导数与泰勒展开。

2. 泰勒展开。 \(f(x+h)=f(x)+f'(x)h+\tfrac12f''(x)h^2+\cdots\)。它是本章分析误差的唯一工具：差分的误差就是被丢掉的那些高阶项。你熟悉的"价格变动 ≈ −久期×Δy + ½×凸性×Δy²"就是债券价格对收益率的二阶泰勒展开。多元版本里 \(h\) 换成向量 \(p\)，一阶项变成 \(\nabla f(x)^Tp\)，二阶项变成 \(\tfrac12p^T\nabla^2f(x)p\)。见 第 00 册第 02 章 和 第 05 章 多元微积分与优化。

3. 链式法则。 若 \(z=g(y)\)、\(y=h(x)\)，则 \(\frac{dz}{dx}=\frac{dz}{dy}\cdot\frac{dy}{dx}\)。多个中间变量时，各条路径的乘积相加。例：\(z=e^{2x}\)，令 \(y=2x\)，则 \(dz/dx=e^y\cdot2=2e^{2x}\)。自动微分就是把链式法则机械地、一步一步地执行。见 第 00 册第 02 章。

4. 梯度、Jacobian、Hessian。 梯度 \(\nabla f\) 是标量函数对每个变量的偏导排成的向量（相当于"所有 Greeks 排成一列"）；Jacobian \(J\) 是向量函数的每个分量对每个变量求偏导排成的矩阵（第 \(j\) 行是第 \(j\) 个输出的梯度）；Hessian \(\nabla^2f\) 是二阶偏导组成的对称矩阵（相当于 Gamma、Vanna、Volga 等二阶 Greeks 的整张表）。\(e_i\) 表示第 \(i\) 个分量为 1、其余为 0 的单位向量，\(x+\epsilon e_i\) 就是"只把第 \(i\) 个变量动 \(\epsilon\)"。见 第 00 册第 05 章。

5. 大 O 记号。 \(O(\epsilon^2)\) 读作"与 \(\epsilon^2\) 同量级"，意思是这一项不超过某个常数乘以 \(\epsilon^2\)。\(\epsilon\) 缩小 10 倍，\(O(\epsilon)\) 的误差缩小 10 倍，\(O(\epsilon^2)\) 的误差缩小 100 倍。见 第 00 册第 07 章 概率中的分析工具。

怎么读这一章

核心必读是 7.2.1（步长与误差，直接决定你算 Greeks 时步长怎么取）、7.3.1–7.3.3（计算图、前向模式、反向模式，AAD 的原理）和 7.4 的实战与解读。7.2.2 与 7.2.4 的图着色第一次读只需抓住"互不干扰的列可以一起扰动"这一个想法，着色规则的细节可以跳过。7.3.4–7.3.6 讲 Jacobian 和 Hessian 的 AD 计算，第一次读看结论（各种量分别是什么代价）即可，等读到第 10、11 章需要 Jacobian 时再回来。7.3.7 的两个局限很短，务必读，它们是实务里真实踩过的坑。

建议顺序：7.1 → 7.2.1 → 7.4 Part A、B → 7.3.1–7.3.3 → 7.4 Part C → 其余各节。

7.1 引言：导数从哪里来

多数非线性优化算法都需要导数：线搜索要梯度，牛顿法要 Hessian 或 Hessian–向量积，非线性最小二乘要 Jacobian。函数简单时可以手工推导，让用户提供；函数复杂时（比如一个几百行的定价程序），就需要自动计算或近似。三类方法：

• 有限差分（finite differencing）：源于 Taylor 定理，看函数在小扰动下的变化来估计导数，例如中心差分 \(\dfrac{\partial f}{\partial x_i}\approx\dfrac{f(x+\epsilon e_i)-f(x-\epsilon e_i)}{2\epsilon}\)。
• 自动微分（automatic differentiation, AD）：把计算函数的程序分解为一系列基本运算的复合，对每一步应用链式法则。有的工具把函数代码变换成同时计算函数与导数的新代码（原书提到的 ADIFOR），有的在给定点运行时记录基本运算、再处理记录得到导数（ADOL-C）。今天的 PyTorch、JAX、TensorFlow 都属于后者的思路。
• 符号微分（symbolic differentiation）：用 Mathematica、Maple 等对代数表达式做符号运算。

本章讨论前两种。导数还有另一个重要用途：最优后的敏感性分析（post-optimal sensitivity analysis），即最优解对参数或约束值小扰动的敏感程度——对量化来说，"组合权重对预期收益输入有多敏感"往往比最优权重本身更重要。

7.2 有限差分

很多优化软件在用户不提供导数时会自动做有限差分。它给出的是近似值，但在很多情形下已经足够。

7.2.1 近似梯度：截断误差与舍入误差

前向差分（forward difference，又称单侧差分）：

整个梯度需要 \(n+1\) 次函数求值。依据是 Taylor 定理：

若 \(\|\nabla^2f\|\le L\)，则 \(|f(x+p)-f(x)-\nabla f(x)^Tp|\le(L/2)\|p\|^2\)（7.3）；取 \(p=\epsilon e_i\)：

推导拆解：从 (7.2) 到 (7.4) 只有三步。

1. (7.2) 是带"拉格朗日余项"的泰勒展开：二阶项里的 Hessian 不在 \(x\) 处取值，而是在 \(x\) 与 \(x+p\) 之间某一点 \(x+tp\)（\(0<t<1\)）取值，这样等式是精确的，不需要写"\(\approx\)"。
2. \(\|\nabla^2f\|\le L\) 的意思是 Hessian 的"大小"（矩阵范数）不超过 \(L\)，于是 \(|\tfrac12p^T\nabla^2f(\cdot)p|\le\tfrac L2\|p\|^2\)，得到 (7.3)。
3. 取 \(p=\epsilon e_i\)：\(\nabla f(x)^Tp=\epsilon\,\partial f/\partial x_i\)（只剩第 \(i\) 个分量），\(\|p\|^2=\epsilon^2\)。把 (7.3) 两边同除以 \(\epsilon\)，就得到"差分商与真导数之差不超过 \(\tfrac L2\epsilon\)"，即 (7.4)。 直观地说，\(L\) 衡量函数的弯曲程度。函数越弯（Gamma 越大），用一条割线去近似切线的误差越大；步长越小，割线越接近切线。

截断误差随 \(\epsilon\) 减小而减小，但计算机有舍入误差。双精度的单位舍入 \(u\approx1.1\times10^{-16}\)（NumPy 中 np.finfo(float).eps 给出 \(2.2\times10^{-16}\)，是 \(2u\)）是每次浮点运算相对误差的上界。粗略地假设计算出的 \(f\) 值相对误差不超过 \(u\)：\(|\mathrm{comp}(f(x))-f(x)|\le uL_f\)（\(L_f\) 是 \(|f|\) 的界）。那么差分商的总误差界为

第一项随 \(\epsilon\) 减小，第二项随 \(\epsilon\) 增大，在 \(\epsilon^2=4L_fu/L\) 时取最小。问题尺度良好时 \(L_f/L\) 适中，取

近乎最优（很多软件采用这个值，实际中常再乘以 \(\max(|x_i|,1)\) 以适应变量的量级），总误差约为 \(\sqrt u\approx10^{-8}\)——前向差分只能得到大约一半的有效数字。

推导拆解：舍入误差项和最优步长是这样来的。

• 计算机存一个数只保留约 16 位有效数字，单位舍入 \(u\) 就是"相对误差的上限"。算出的 \(f(x+\epsilon e_i)\) 和 \(f(x)\) 各自可能偏差 \(uL_f\)，两者相减最坏偏差 \(2uL_f\)，再除以 \(\epsilon\)，就是 (7.5) 的第二项 \(2uL_f/\epsilon\)。问题的根源是两个几乎相等的数相减：例如 \(100.000000123-100.000000000\)，前面 9 位全部抵消，结果只剩后面几位可信，而后面几位恰恰是舍入误差所在的位置。
• 求 (7.5) 的最小值：对 \(\epsilon\) 求导并令其为零，\(\tfrac L2-\tfrac{2uL_f}{\epsilon^2}=0\)，得 \(\epsilon^2=4uL_f/L\)。代回去，两项相等，总误差 \(=2\sqrt{uL_fL}\)。当 \(L_f\)、\(L\) 都在 1 的量级时，\(\epsilon^*\approx2\sqrt u\)、误差 \(\approx2\sqrt u\)，所以说"\(\epsilon=\sqrt u\) 近乎最优、精度约 \(10^{-8}\)"。
• 这和"偏差–方差权衡"同构：截断误差像偏差（步长越大越偏），舍入误差像方差（步长越小，噪声被 \(1/\epsilon\) 放大得越厉害），最优点在两者相等处。

金融直觉：风险系统里 DV01 通常用 1bp（\(10^{-4}\)）的扰动，而不是 \(10^{-8}\)。原因有两个：一是收益率本身就在 \(10^{-2}\) 量级，1bp 已经是相对很小的扰动；二是定价程序内部常有插值、迭代求解等"噪声"，实际的 \(L_f\) 误差远大于 \(u\)，按 (7.5) 的逻辑，噪声越大，最优步长越大。本章 7.4 的 Part B 会把这一点说透。

中心差分（central difference）：

需要 \(2n\) 次函数求值，代价约为前向差分的两倍。若 Hessian Lipschitz 连续，

取 \(p=\pm\epsilon e_i\) 相减，二阶项抵消，截断误差为 \(O(\epsilon^2)\)。把舍入误差一起考虑，最优步长约为 \(\epsilon\approx u^{1/3}\)，可达精度约 \(u^{2/3}\approx10^{-11}\)（原书习题 7.1；PDF 正文此处印作"\(\epsilon=u^{2/3}\)"，与习题不一致，以习题为准）。多出的几位有效数字有时值得额外的代价。

推导拆解：为什么二阶项会抵消。写出两边展开（一元记号，\(g=\partial f/\partial x_i\)，\(H=\partial^2f/\partial x_i^2\)）： \(f(x+\epsilon e_i)=f+g\epsilon+\tfrac12H\epsilon^2+\tfrac16T\epsilon^3+\cdots\)， \(f(x-\epsilon e_i)=f-g\epsilon+\tfrac12H\epsilon^2-\tfrac16T\epsilon^3+\cdots\)。 两式相减：\(f\) 和 \(\tfrac12H\epsilon^2\) 都抵消，剩 \(2g\epsilon+\tfrac13T\epsilon^3\)。除以 \(2\epsilon\) 得 \(g+\tfrac16T\epsilon^2\)，所以误差是 \(O(\epsilon^2)\)。 最优步长：总误差 \(\approx\tfrac M6\epsilon^2+\tfrac{uL_f}{\epsilon}\)，求导令其为零得 \(\epsilon^3\propto u\)，即 \(\epsilon\propto u^{1/3}\approx6\times10^{-6}\)；代回得误差 \(\propto u^{2/3}\approx10^{-11}\)。 金融上，这就是为什么算 Delta 时"上下各 bump 一次"比"只往上 bump"准得多：单边 bump 的误差里混着 Gamma 的贡献，对称 bump 把 Gamma 的贡献对消了。

7.2.2 近似稀疏 Jacobian

对向量函数 \(r:\mathbb{R}^n\to\mathbb{R}^m\)（非线性最小二乘的残差、非线性方程组），Jacobian \(J(x)\) 的第 \(j\) 行是 \(\nabla r_j(x)^T\)。由 Taylor 定理，\(\|r(x+p)-r(x)-J(x)p\|\le(L/2)\|p\|^2\)（7.9）。

• Jacobian–向量积（例如非线性方程组的非精确牛顿法所需）只要一次额外求值：
  \[J(x)p\approx\frac{r(x+\epsilon p)-r(x)}{\epsilon},\tag{7.10}\]
  精度 \(O(\epsilon)\)；也可以用双侧版本。
• 完整 Jacobian 逐列计算：
  \[\frac{\partial r}{\partial x_i}(x)\approx\frac{r(x+\epsilon e_i)-r(x)}{\epsilon},\tag{7.11}\]
  需要 \(n+1\) 次求值。但若 \(J\) 稀疏，可以同时估计多列，有时只需三四次。

例（原书 (7.12)）

每个分量只依赖两三个相邻变量，Jacobian 是三对角的。扰动 \(\epsilon e_1\) 只影响 \(r_1,r_2\)；扰动 \(\epsilon e_4\) 只影响 \(r_3,r_4,r_5\)——两者互不干扰。所以取 \(p=\epsilon(e_1+e_4)\) 做一次求值，就能同时得到第 1 列的 \((1,1),(2,1)\) 元和第 4 列的 \((3,4),(4,4),(5,4)\) 元：

（PDF 文本中第二组的下标印作 \(\partial r_4/\partial x_3\) 等，并误称为"Hessian 第四列"，按上下文应为 Jacobian 第 4 列，此处已更正。）同样用 \(\epsilon(e_2+e_5)\)、\(\epsilon(e_3+e_6)\)，总共 3 次额外求值就得到整个 Jacobian。对任意 \(n\) 都只需 3 个扰动向量：\(\epsilon(e_1+e_4+e_7+\cdots)\)、\(\epsilon(e_2+e_5+\cdots)\)、\(\epsilon(e_3+e_6+\cdots)\)——同组的列在任何一行都不会同时非零。

图着色（graph coloring）把这一想法一般化。构造列关联图（column intersection graph）\(\mathcal{G}\)：\(n\) 个节点对应 \(n\) 列，若某个 \(r_j\) 同时依赖 \(x_i\) 和 \(x_k\)（即第 \(i,k\) 列在某一行都非零）就在 \(i,k\) 之间连边（原书图 7.1）。给节点着色，使相邻节点颜色不同；每种颜色对应一个扰动向量 \(p=\epsilon(e_{i_1}+\cdots+e_{i_\ell})\)。求最少着色数是困难问题，但有便宜的近似最优算法（Curtis–Powell–Reid，Coleman–Moré）。用更一般的扰动向量，额外求值次数可以不超过每行非零元的最大个数（Newsam–Ramsdell）。

白话解释：着色的道理是"互不相干的扰动可以同时做"。如果第 1 列和第 4 列在任何一行都不会同时非零，那就说明没有一个输出同时依赖 \(x_1\) 和 \(x_4\)；同时扰动 \(x_1\) 和 \(x_4\) 后，每个输出的变化只能来自其中一个，不会混在一起，各自读出即可。"颜色"只是给"可以一起扰动的那一组变量"起的名字，颜色数就是需要额外求值的次数。 金融上对应的场景：曲线自举时，3 个月存款只依赖曲线最前端的节点，10 年互换主要依赖中长端节点。如果两个节点不出现在同一个工具的定价里，就可以把它们一起 bump，用一次重新定价同时读出两组敏感度。

7.2.3 近似 Hessian

• 有梯度、没有 Hessian：由 \(\nabla f(x+p)=\nabla f(x)+\nabla^2f(x)p+O(\|p\|^2)\)（7.18），
  \[\nabla^2f(x)e_i\approx\frac{\nabla f(x+\epsilon e_i)-\nabla f(x)}{\epsilon},\tag{7.19}\]
  误差 \(O(\epsilon)\)，完整 Hessian 需要 \(n+1\) 次梯度求值。逐列得到的近似矩阵不一定对称，可取 \((H+H^T)/2\)。
• Hessian–向量积（Newton–CG 所需）只要一次额外梯度：
  \[\nabla^2f(x)p\approx\frac{\nabla f(x+\epsilon p)-\nabla f(x)}{\epsilon};\tag{7.20}\]
  再算一次 \(\nabla f(x-\epsilon p)\) 就得到中心差分版本。
• 连梯度都没有：用 (7.8)，取 \(p=\epsilon e_i,\ \epsilon e_j,\ \epsilon(e_i+e_j)\) 组合：
  \[\frac{\partial^2f}{\partial x_i\partial x_j}(x)=\frac{f(x+\epsilon e_i+\epsilon e_j)-f(x+\epsilon e_i)-f(x+\epsilon e_j)+f(x)}{\epsilon^2}+O(\epsilon),\tag{7.21}\]
  完整 Hessian 需要 \(n(n+1)/2+n\) 个点的函数值；稀疏时跳过已知为零的元素。注意分母是 \(\epsilon^2\)，舍入误差被放大得更厉害，最优步长约为 \(u^{1/3}\) 而非 \(\sqrt u\)。

7.2.4 近似稀疏 Hessian

Hessian 是 \(\nabla f\) 的 Jacobian，可以直接套用稀疏 Jacobian 技术，但那样忽略了对称性：估计出 \((i,j)\) 元也就有了 \((j,i)\) 元，利用这一点可以大幅节省。

例（原书 (7.22)）

它的 Hessian 是"箭头形"的：第一行、第一列和对角线非零。由于第一行处处非零，列关联图是完全图，按 Jacobian 规则需要 \(n+1\) 次梯度求值。利用对称性：先用 \(p=\epsilon e_1\) 估计第一列（也就得到了第一行）；剩下的只有对角元 \(2,\dots,n\)，它们之间互不相连，可以着同一种颜色：

\(\nabla f\) 的第 \(i\) 个分量（\(i\ge2\)）只受 \(x_i\)（以及 \(x_1\)，但 \(x_1\) 未被扰动）的影响，所以 \(\dfrac{\partial^2f}{\partial x_i^2}\approx\dfrac{[\nabla f(x+p)-\nabla f(x)]_i}{\epsilon}\)。总共只需在 \(x\) 和另外两点算梯度。

一般方法用邻接图（adjacency graph：\(i\ne k\) 且 \(\partial^2f/\partial x_i\partial x_k\ne0\) 时连边），着色要求：相邻节点颜色不同，且任何长度为 3 的路径（\(i_1-i_2-i_3-i_4\)）至少用三种颜色（Coleman–Moré）。

7.3 自动微分

自动微分利用函数的计算过程得到导数的精确值（在浮点精度内，不是近似）。它的基础是两点：任何函数的计算都可以分解为一列一元或二元基本运算（加、减、乘、除、幂，三角、指数、对数）；而链式法则

告诉我们如何把基本运算的导数组合起来。两种基本模式：前向模式和反向模式。

7.3.1 计算图

例（原书 (7.26)）

引入中间变量：

把每个变量画成节点，若 \(x_j\) 的计算直接用到 \(x_i\) 就画一条有向边 \(i\to j\)，得到计算图（原书图 7.2）。称 \(i\) 为 \(j\) 的父节点、\(j\) 为 \(i\) 的子节点。父节点的值都已知时就能算出子节点，计算从左向右流动，称为前向扫描（forward sweep）。AD 工具会自动识别中间量、构造计算图，用户无需手工分解。注意 \(x_4\) 被 \(x_6\) 和 \(x_7\) 共用——这种公共子表达式是 AD 高效的来源之一。

7.3.2 前向模式

前向模式在计算每个中间变量的同时，计算它沿某个方向 \(p\) 的方向导数：

目标是 \(D_px_9=\nabla f(x)^Tp\)。在自变量处 \(D_px_i=p_i\)，\(p\) 称为种子向量（seed vector）。每个中间变量的方向导数由其父节点按链式法则得到，例如 \(x_7=x_4x_5\)：

实现上，软件给每个标量 \(w\) 附带一个 \(D_pw\)，每做一次运算就同步做相应的链式运算，例如 \(z=w/y\)：

白话解释：前向模式就是"每个数都随身带一个小账本"。普通计算只记 \(w\) 的值；前向模式同时记"如果输入沿方向 \(p\) 动一点，\(w\) 会跟着动多少"，即 \(D_pw\)。每做一步运算，值按原规则算，账本按链式法则更新。(7.30) 就是商的求导法则 \((w/y)'=w'/y-wy'/y^2\)，只是把"\('\)"换成了"沿 \(p\) 方向的变化率"。 例：取 \(p=e_1\)（只动 \(x_1\)），在 \(x=(1,2,\pi/2)\) 处，\(D_px_1=1\)、\(D_px_2=0\)；\(x_4=x_1x_2\) 的账本为 \(x_2D_px_1+x_1D_px_2=2\)，表示"\(x_1\) 动 1，\(x_4\) 动 2"。一路传到 \(x_9\)，账本里就是 \(\partial f/\partial x_1\)。 一次前向传播只能得到"一个方向"的导数，相当于一次 bump；要得到全部 \(n\) 个偏导，就得做 \(n\) 次（或带 \(n\) 个账本），所以代价与有限差分同量级，好处只是没有步长问题、结果精确。

不必同时存储所有节点的值：一个节点的所有子节点算完后，它就可以被覆盖。实现方式有两种：预编译器（把源代码变换为扩展代码），或在 C++ 等语言中运算符重载（今天的 Python AD 库也多用此法）。

要得到完整梯度，就对 \(n\) 个种子 \(e_1,\dots,e_n\) 同时做前向传播。代价可能显著：一次除法要引起约 \(2n\) 次乘法和 \(n\) 次加法，存储也可能增加 \(n\) 倍。计算早期很多 \(D_{e_j}x_i\) 为零，可以用稀疏数据结构节省。前向模式求完整梯度的代价大约是函数求值的 \(n\) 倍——与有限差分同一量级，但结果是精确的。

7.3.3 反向模式

反向模式先完成函数的计算（前向扫描），然后反向扫描计算图，求出 \(f\) 对每个变量（自变量和中间变量）的偏导数，最后在自变量节点处组装出梯度。

• 每个节点关联一个标量伴随变量（adjoint variable）\(\bar x_i\)，初始化为零；最右边的输出节点置 \(\bar x_N=1\)（因为 \(\partial f/\partial x_N=1\)）。
• 链式法则写成
  \[\frac{\partial f}{\partial x_i}=\sum_{j\text{ 为 }i\text{ 的子节点}}\frac{\partial f}{\partial x_j}\frac{\partial x_j}{\partial x_i},\tag{7.31}\]
  每当某一项已知就累加进去：
  \[\bar x_i\mathrel{+}=\frac{\partial f}{\partial x_j}\frac{\partial x_j}{\partial x_i}.\tag{7.32}\]
  一个节点的所有子节点都贡献完毕后，\(\bar x_i=\partial f/\partial x_i\)，称该节点"完成"（finalized），然后它再向自己的父节点贡献。计算方向是从子到父，与求值方向相反。
• 反向扫描用的是数值而不是公式：前向扫描时不仅计算 \(x_i\)，还要计算并存储每条边上的偏导数值 \(\partial x_j/\partial x_i\)。

白话解释：前向模式问的是"这个输入动一下，会影响哪些下游量"；反向模式反过来问"最终输出对每个中间量有多敏感"。伴随变量 \(\bar x_i=\partial f/\partial x_i\) 就是"中间量 \(x_i\) 每变 1 单位，最终结果变多少"。 (7.31) 说的是：\(x_i\) 只能通过它的子节点影响 \(f\)，所以 \(x_i\) 的敏感度 = 各子节点的敏感度 × 该子节点对 \(x_i\) 的局部导数，再求和。因为要用到子节点的敏感度，就必须从输出端往回算。 金融直觉：这很像损益归因倒着做。总损益对一个衍生品账簿的敏感度已知（\(\bar x_N=1\)），先算它对每个交易的敏感度，再算每个交易对各个风险因子的敏感度，层层往下分摊。每个风险因子最后收到的，就是来自所有路径的贡献之和。

数值例（原书图 7.3）：\(x=(1,2,\pi/2)^T\)。前向扫描得 \(x_4=2\)，\(x_5=1\)，\(x_6=e^2\)，\(x_7=2\)，\(x_8=2+e^2\)，\(x_9=(4+2e^2)/\pi\)；边上的偏导数值：\(\frac{\partial x_4}{\partial x_1}=2\)，\(\frac{\partial x_4}{\partial x_2}=1\)，\(\frac{\partial x_5}{\partial x_3}=\cos\frac\pi2=0\)，\(\frac{\partial x_6}{\partial x_4}=e^2\)，\(\frac{\partial x_7}{\partial x_4}=x_5=1\)，\(\frac{\partial x_7}{\partial x_5}=x_4=2\)，\(\frac{\partial x_8}{\partial x_6}=\frac{\partial x_8}{\partial x_7}=1\)，\(\frac{\partial x_9}{\partial x_8}=\frac1{x_3}=\frac2\pi\)，\(\frac{\partial x_9}{\partial x_3}=-\frac{x_8}{x_3^2}=-\frac{2+e^2}{(\pi/2)^2}\)。

反向扫描：\(\bar x_9=1\)；节点 9 向父节点 3 和 8 贡献

节点 3 还在等待子节点 5 的贡献；节点 8 只有子节点 9，已完成，于是 \(\bar x_6\mathrel{+}=\frac2\pi\)，\(\bar x_7\mathrel{+}=\frac2\pi\)，依此继续更新节点 4、5……最终

（验证第一个分量：\(\partial f/\partial x_1=x_2(\sin x_3+e^{x_1x_2})/x_3=2(1+e^2)/(\pi/2)=(4+4e^2)/\pi\)。节点 5 对 \(\bar x_3\) 的贡献是 \(\bar x_5\cdot\cos(\pi/2)=0\)，所以 \(\bar x_3\) 就是 (7.33) 的值。）

反向模式的关键优点：对标量函数 \(f:\mathbb{R}^n\to\mathbb{R}\)，计算梯度的额外算术量至多为函数求值的 4–5 倍，与 \(n\) 无关。原书的论证：每个基本运算在反向扫描中只引起常数次运算（例如 (7.33) 中节点 9 的处理需要两次乘法、一次除法、一次加法）。前向模式则可能需要 \(n\) 倍。对向量函数 \(r:\mathbb{R}^n\to\mathbb{R}^m\)，随着 \(m\) 增大，两种模式的代价趋于接近。

金融直觉：选哪种模式，看"输入多还是输出多"。前向模式一次处理一个输入方向，反向模式一次处理一个输出。一个奇异期权的价格（1 个输出）依赖几百个曲线节点和波动率格点（几百个输入），反向模式一次扫描就拿到全部敏感度；有限差分或前向模式则要重新定价几百次。反过来，如果是"1 个参数影响 1000 个产品的价格"，就用前向模式。 量级感受：若一次定价 0.1 秒、输入 500 个，中心差分要 \(2\times500\times0.1=100\) 秒；AAD 约 \(0.1\times5=0.5\) 秒。这就是 2010 年以后各大行在 XVA 和监管资本计算中大规模上 AAD 的原因。

缺点是存储：需要保存整个计算图（每个基本运算一个节点：中间结果、指向一两个父节点的指针、边上的偏导数）。反向扫描时按写入的逆序读取，访问模式简单，可以通过运算符重载实现（ADOL-C 的做法），反向扫描只是一次函数调用。但存储可能巨大：原书估计每节点 20 字节，在一台每秒一亿次浮点运算的机器上求值 1 秒的函数，计算图就可达 2 GB。解决办法是检查点（checkpointing）：只保存图的若干片段的起点，反向时对每一片段重新做局部的前向扫描，以额外计算换存储。

7.3.4 向量函数与部分可分性

非线性最小二乘与非线性方程中 \(r:\mathbb{R}^n\to\mathbb{R}^m\)，计算图最右边有 \(m\) 个没有子节点的输出节点，Jacobian 为

部分可分（partially separable）函数是指

其中每个元素函数（element function）\(f_i\) 只依赖 \(x\) 的少数分量。令 \(r(x)=(f_1(x),\dots,f_{n_e}(x))^T\)，则

\(J\) 的多数列非零元很少，可以用图着色高效计算再恢复梯度。很多大规模问题天然是部分可分的（本册第 9 章有对应的拟牛顿方法）。约束优化中，把目标和约束放在一起计算（\(r(x)=(f(x),c_1(x),\dots)\)）可以共享公共子表达式，减少总工作量。

7.3.5 用自动微分计算 Jacobian

• 前向模式：给种子 \(p\)，输出节点得到 \(D_pr_j=(\nabla r_j)^Tp\)，合起来就是 Jacobian–向量积 \(J(x)p\)。完整 Jacobian 取 \(p=e_1,\dots,e_n\)；稀疏时用着色选种子。算术量约为函数的"种子数"倍。
• 反向模式：选种子 \(q\in\mathbb{R}^m\)，对标量函数 \(r(x)^Tq\) 做反向扫描，得到 \(\nabla[r(x)^Tq]=J(x)^Tq\)，即 Jacobian 转置–向量积。完整 Jacobian 取 \(q=e_1,\dots,e_m\)；稀疏时对 \(J^T\) 着色。算术量不超过种子数的 5 倍，存储不超过标量情形。
• 两种模式可以组合：前向种子得到一部分列，反向种子得到其余的行。
• 非精确牛顿法等只需反复计算 \(J(x)p\)，一次前向扫描即可，代价与一次函数求值相当。

7.3.6 用自动微分计算 Hessian

前向模式：对种子对 \((p,q)\) 定义二阶量

在前向扫描中与 \(x_i\)、\(D_px_i\)、\(D_qx_i\) 同步计算；自变量处 \(D_{pq}=0\)；输出节点得到 \(p^T\nabla^2f(x)q\)。传播规则：加法 \(x_i=x_j+x_k\) 时

一元运算 \(x_i=L(x_j)\) 时

推导拆解：(7.39) 的二阶规则就是"对链式法则再求一次导"。先看一元情形 \(x_i=L(x_j(x))\)：

1. 一阶：\(\nabla x_i=L'(x_j)\nabla x_j\)，两边与 \(p\) 做内积就是 \(D_px_i=L'D_px_j\)。
2. 再对 \(x\) 求导（乘积法则）：\(\nabla^2x_i=L''(x_j)\nabla x_j\nabla x_j^T+L'(x_j)\nabla^2x_j\)。第一项来自"\(L'(x_j)\) 本身也随 \(x\) 变"，用了链式法则；第二项来自"\(\nabla x_j\) 随 \(x\) 变"。
3. 左乘 \(p^T\)、右乘 \(q\)：\(p^T\nabla x_j=D_px_j\)，\(\nabla x_j^Tq=D_qx_j\)，\(p^T\nabla^2x_jq=D_{pq}x_j\)，即得 (7.39)。 一元版本你很熟：\(\frac{d^2}{dx^2}L(g(x))=L''(g)g'^2+L'(g)g''\)。例如债券价格 \(P=e^{-yT}\) 中 \(y\) 又是某个参数的函数时，凸性项就分成"指数函数本身的弯曲"和"\(y\) 自身的弯曲"两部分。

• 稠密 Hessian：对所有单位向量对 \((e_j,e_k)\)，\(k\le j\)，共 \(n(n+1)/2\) 对（由对称性只算下三角）；已知稀疏结构时只算可能非零的位置。
• 只要 Hessian–向量积 \(\nabla^2f(x)q\)：计算 \(D_{e_j}x_i\)、\(D_qx_i\) 和 \(D_{e_jq}x_i\)（\(j=1..n\)），输出节点给出 \(e_j^T\nabla^2f(x)q\)；代价约为 \(2n\) 的小倍数。
• 另一种更简单的做法只用"单方向"的一、二阶导：
  \[[\nabla^2f(x)]_{ij}=\tfrac12\left[(e_i+e_j)^T\nabla^2f(e_i+e_j)-e_i^T\nabla^2fe_i-e_j^T\nabla^2fe_j\right],\tag{7.40}\]
  只需对 \(p=e_i,e_j,e_i+e_j\) 传播 \(D_p\)、\(D_{pp}\)，不需要交叉项 \(D_{pq}\)。\(D_pf\)、\(D_{pp}f\) 就是一元函数 \(\psi(t)=f(x+tp)\)（7.41）在 \(t=0\) 处的一、二阶导数，这种思路可以推广到更高阶导数。

反向模式求 Hessian–向量积：先用前向模式同时算出 \(f\) 与 \(\nabla f(x)^Tq\)（累积 \(x_i\)、\(D_qx_i\)），再对"计算 \(\nabla f(x)^Tq\) 的程序"做一次标准的反向扫描，自变量节点处得到

算术量与 \(n\) 无关：前向部分是小倍数，反向再乘至多 5，原书估计总共约为 \(f\) 的 12 倍。完整 Hessian 取 \(q=e_1,\dots,e_n\)，至多 \(12n\) 倍；有稀疏结构时用着色，约为 \(12N_c\) 倍（\(N_c\) 为种子数）。这正是第 6 章 Hessian-free Newton–CG 所需的工具，也是现代深度学习框架中"Hessian–向量积"的标准实现（前向套反向）。

7.3.7 局限

• 截断误差的导数：若 \(f\) 的计算本身依赖于数值近似（如 PDE 的数值解），计算值是 \(\hat f(x)=f(x)+\tau(x)\)。\(|\tau|\) 很小，但 \(\tau'\) 未必小，所以 AD 给出的"精确导数"可能与真实导数相差很大（有限差分同样受影响）。用分段有理函数近似三角函数时也有类似问题。

金融直觉：二叉树定价是最典型的例子。\(N\) 步二叉树的价格 \(\hat f(S)\) 与 Black–Scholes 价格之差 \(\tau(S)\) 很小，但随 \(S\) 呈锯齿状来回振荡（节点相对执行价的位置在变）。锯齿的幅度小，斜率却不小，所以对树价格直接求导（无论 AD 还是差分）得到的 Delta、尤其是 Gamma，会出现明显的抖动。实务中常用的办法是对树做平滑（如在倒数第二步用 BS 公式代替），或直接从树的节点上读 Delta，而不是对整个定价函数求导。 同样，带数字期权、障碍期权这类不连续 payoff 的蒙特卡洛，对每条路径做 AD 得到的导数几乎处处为零（payoff 是阶梯函数），也属于这一类问题，需要先光滑化。

• 代码中的分支：病态例子——\(f(x)=x-1\) 被写成

if (x == 1.0) then f = 0.0 else f = x − 1.0

AD 在 \(x=1\) 处只看到了常数分支，给出 \(f'(1)=0\)，而正确值是 1。

结论：AD 是日益成熟的技术，让优化算法能用于复杂函数，也便于解释最优解（敏感性分析）；但它不是万灵药，用户仍需思考导数是怎么算的。

7.4 量化实战：Greeks 的差分、公共随机数与 AAD

期权 Greeks 是"价格对市场参数的导数"，是定价与对冲中最常计算的导数。本节分四部分：

• A：用 Black–Scholes 看涨期权的 Delta 验证 (7.5)–(7.7) 的步长理论。
• B：蒙特卡洛定价时，函数值带统计噪声（远大于 \(u\)），差分必须使用公共随机数（common random numbers）。
• C：用几十行代码实现一个反向模式 AD（运算符重载 + 计算图 + 伴随累加），先复现原书例 (7.26)，再用一次反向扫描得到 Black–Scholes 价格对 \(S,\sigma,r,T\) 的全部一阶敏感度——这就是衍生品行业所说的伴随算法微分（Adjoint Algorithmic Differentiation, AAD）的原理；最后演示分支陷阱。
• D：用 3 种颜色的扰动向量，以 3 次额外求值得到原书例 (7.12) 的三对角 Jacobian。

import numpy as np
from math import erf, sqrt, pi
from scipy.stats import norm

# ---------- Part A: 有限差分步长与误差（Black–Scholes 看涨期权 Delta） ----------
def bs_call(S, K=100.0, T=0.5, r=0.03, sigma=0.25):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d1 - sigma * np.sqrt(T))
S0 = 100.0
d1 = (np.log(S0 / 100) + (0.03 + 0.5 * 0.25**2) * 0.5) / (0.25 * np.sqrt(0.5))
delta_true = norm.cdf(d1)
print("相对步长 ε    前向差分误差   中心差分误差")
for e in [1e-2, 1e-4, 1e-5, 1e-6, 1e-8, 1e-10, 1e-12]:
    h = e * S0
    fwd = (bs_call(S0 + h) - bs_call(S0)) / h
    ctr = (bs_call(S0 + h) - bs_call(S0 - h)) / (2 * h)
    print("  %.0e      %.2e       %.2e" % (e, abs(fwd - delta_true), abs(ctr - delta_true)))
u = np.finfo(float).eps
print("理论最优：前向 ε≈√u=%.1e，中心 ε≈u^(1/3)=%.1e" % (np.sqrt(u), u ** (1 / 3)))

# ---------- Part B: 蒙特卡洛 Greeks —— 噪声下的差分，必须用公共随机数 ----------
def mc_call(S, Z, K=100.0, T=0.5, r=0.03, sigma=0.25):
    ST = S * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
    return np.exp(-r * T) * np.maximum(ST - K, 0).mean()
rng = np.random.default_rng(7); N, h = 100_000, 1.0
est_crn, est_ind = [], []
for _ in range(40):
    Z = rng.standard_normal(N)
    est_crn.append((mc_call(S0 + h, Z) - mc_call(S0 - h, Z)) / (2 * h))              # 公共随机数
    est_ind.append((mc_call(S0 + h, rng.standard_normal(N)) - mc_call(S0 - h, Z)) / (2 * h))  # 独立随机数
print("\nMC Delta（真值 %.4f）：公共随机数 均值 %.4f 标准差 %.4f；独立随机数 均值 %.4f 标准差 %.4f" %
      (delta_true, np.mean(est_crn), np.std(est_crn), np.mean(est_ind), np.std(est_ind)))

# ---------- Part C: 一个最小的反向模式自动微分（运算符重载 + 计算图） ----------
class Var:
    def __init__(self, val, parents=()):
        self.val, self.parents, self.adj = val, parents, 0.0     # parents: [(父节点, ∂self/∂父节点)]
    def __add__(a, b): b = b if isinstance(b, Var) else Var(b); return Var(a.val + b.val, [(a, 1.0), (b, 1.0)])
    def __sub__(a, b): b = b if isinstance(b, Var) else Var(b); return Var(a.val - b.val, [(a, 1.0), (b, -1.0)])
    def __mul__(a, b): b = b if isinstance(b, Var) else Var(b); return Var(a.val * b.val, [(a, b.val), (b, a.val)])
    def __truediv__(a, b):
        b = b if isinstance(b, Var) else Var(b)
        return Var(a.val / b.val, [(a, 1 / b.val), (b, -a.val / b.val**2)])
    __radd__ = __add__; __rmul__ = __mul__
    def __rsub__(a, b): return Var(b) - a
    def __neg__(a): return Var(-a.val, [(a, -1.0)])
def exp(a): v = np.exp(a.val); return Var(v, [(a, v)])
def log(a): return Var(np.log(a.val), [(a, 1 / a.val)])
def sin(a): return Var(np.sin(a.val), [(a, np.cos(a.val))])
def vsqrt(a): v = sqrt(a.val); return Var(v, [(a, 0.5 / v)])
def Phi(a): return Var(0.5 * (1 + erf(a.val / sqrt(2))), [(a, np.exp(-0.5 * a.val**2) / sqrt(2 * pi))])

def backward(out):
    """反向扫描：按拓扑序的逆序累加伴随变量 x̄_i += x̄_j · ∂x_j/∂x_i  —— (7.32)"""
    order, seen = [], set()
    def topo(v):
        if id(v) in seen: return
        seen.add(id(v))
        for p, _ in v.parents: topo(p)
        order.append(v)
    topo(out); out.adj = 1.0
    for v in reversed(order):
        for p, dp in v.parents: p.adj += v.adj * dp
    return len(order)

# 原书例 (7.26)：f = (x1 x2 sin x3 + exp(x1 x2)) / x3，在 x = (1, 2, π/2)
x1, x2, x3 = Var(1.0), Var(2.0), Var(pi / 2)
f = (x1 * x2 * sin(x3) + exp(x1 * x2)) / x3
backward(f)
e2 = np.exp(2)
print("\n例 7.26 反向模式梯度:", np.round([x1.adj, x2.adj, x3.adj], 6))
print("原书解析结果        :", np.round([(4 + 4 * e2) / pi, (2 + 2 * e2) / pi, (-8 - 4 * e2) / pi**2], 6))

# 一次反向扫描得到 Black–Scholes 价格对 (S, σ, r, T) 的全部一阶敏感度（AAD 的原理）
S, sig, r, T, K = Var(100.0), Var(0.25), Var(0.03), Var(0.5), 100.0
sqT = vsqrt(T)
d1 = (log(S / K) + (r + 0.5 * sig * sig) * T) / (sig * sqT)
price = S * Phi(d1) - K * exp(-r * T) * Phi(d1 - sig * sqT)
nodes = backward(price)
d1v = d1.val; d2v = d1v - 0.25 * np.sqrt(0.5)
analytic = {"Delta": norm.cdf(d1v), "Vega": 100 * norm.pdf(d1v) * np.sqrt(0.5),
            "Rho": 100 * 0.5 * np.exp(-0.015) * norm.cdf(d2v),
            "dV/dT": 100 * norm.pdf(d1v) * 0.25 / (2 * np.sqrt(0.5)) + 0.03 * 100 * np.exp(-0.015) * norm.cdf(d2v)}
print("BS 价格 %.4f，计算图 %d 个节点" % (price.val, nodes))
for name, v in zip(analytic, [S.adj, sig.adj, r.adj, T.adj]):
    print("  %-6s AD = %.6f   解析 = %.6f" % (name, v, analytic[name]))

# AD 的分支陷阱（原书 7.2 节末）：f(x) = x - 1，但代码在 x == 1 处走了另一个分支
def f_branch(x): return Var(0.0) if x.val == 1.0 else x - 1.0
x = Var(1.0); y = f_branch(x); backward(y)
print("分支陷阱：f(x)=x-1 在 x=1 处 AD 给出 f'(1) = %.1f（正确值 1）" % x.adj)

# ---------- Part D: 稀疏 Jacobian 的图着色差分（原书例 (7.12)，n = 6） ----------
def resid(x):
    n = len(x); r = np.zeros(n)
    r[0] = 2 * (x[1]**3 - x[0]**2)
    for i in range(1, n - 1):
        r[i] = 3 * (x[i]**3 - x[i - 1]**2) + 2 * (x[i + 1]**3 - x[i]**2)
    r[n - 1] = 3 * (x[n - 1]**3 - x[n - 2]**2)
    return r
n = 6; x = np.linspace(0.5, 1.5, n); eps = np.sqrt(u) * 1.0; r0 = resid(x)
J = np.zeros((n, n))
for color in range(3):                       # 三种颜色：{1,4}, {2,5}, {3,6}
    cols = np.arange(color, n, 3)
    p = np.zeros(n); p[cols] = eps
    dr = (resid(x + p) - r0) / eps           # 一次额外求值估计同色的所有列
    for c in cols:
        rows = [i for i in (c - 1, c, c + 1) if 0 <= i < n]   # 第 c 列的非零行（三对角结构）
        J[rows, c] = dr[rows]
J_full = np.column_stack([(resid(x + eps * np.eye(n)[i]) - r0) / eps for i in range(n)])
print("\n着色法：3 次额外求值；逐列法：%d 次；两者最大差 %.1e" % (n, np.abs(J - J_full).max()))
print(np.round(J, 3))

运行输出：

相对步长 ε    前向差分误差   中心差分误差
  1e-02      1.10e-02       7.33e-05
  1e-04      1.11e-04       7.34e-09
  1e-05      1.11e-05       7.81e-11
  1e-06      1.11e-06       5.69e-11
  1e-08      2.72e-09       8.31e-10
  1e-10      2.10e-07       2.10e-07
  1e-12      2.04e-05       2.04e-05
理论最优：前向 ε≈√u=1.5e-08，中心 ε≈u^(1/3)=6.1e-06

MC Delta（真值 0.5688）：公共随机数 均值 0.5687 标准差 0.0020；独立随机数 均值 0.5666 标准差 0.0235

例 7.26 反向模式梯度: [10.681278  5.340639 -3.805241]
原书解析结果        : [10.681278  5.340639 -3.805241]
BS 价格 7.7603，计算图 28 个节点
  Delta  AD = 0.568769   解析 = 0.568769
  Vega   AD = 27.789321   解析 = 27.789321
  Rho    AD = 24.558325   解析 = 24.558325
  dV/dT  AD = 8.420830   解析 = 8.420830
分支陷阱：f(x)=x-1 在 x=1 处 AD 给出 f'(1) = 0.0（正确值 1）

着色法：3 次额外求值；逐列法：6 次；两者最大差 0.0e+00
[[-2.    2.94  0.    0.    0.    0.  ]
 [-3.    1.61  4.86  0.    0.    0.  ]
 [ 0.   -4.2   3.69  7.26  0.    0.  ]
 [ 0.    0.   -5.4   6.49 10.14  0.  ]
 [ 0.    0.    0.   -6.6  10.01 13.5 ]
 [ 0.    0.    0.    0.   -7.8  20.25]]

解读

• A：步长。前向差分误差随 \(\epsilon\) 减小线性下降，在 \(\epsilon=10^{-8}\)（\(\approx\sqrt u\)）附近达到最小 \(\sim3\times10^{-9}\)，再减小步长误差反而增大（舍入误差主导）。中心差分误差在 \(\epsilon=10^{-5}\)–\(10^{-6}\)（\(\approx u^{1/3}\)）附近达到最小 \(\sim6\times10^{-11}\)，比前向差分多出约两位有效数字。\(\epsilon=10^{-12}\) 时两者都只剩约 5 位有效数字。这与 (7.5)–(7.7) 的理论完全吻合。
• B：噪声下的差分。蒙特卡洛价格的统计误差约为 \(10^{-2}\) 量级，远大于机器精度，所以 (7.5) 中的 \(uL_f\) 要换成噪声水平，最优步长也要相应增大（这里用了 \(h=1\)，即 1% 的价格扰动）。更关键的是：用同一组随机数计算上下两个价格，Delta 估计的标准差只有 0.002；用独立随机数，标准差高达 0.024，是前者的 12 倍——两次价格估计的噪声不再抵消，而是被 \(1/(2h)\) 放大。实务中用有限差分算蒙特卡洛 Greeks，必须固定随机数种子；对数字期权、障碍期权这类不连续的 payoff，还需要光滑化或似然比方法。
• C：AAD。手写的反向模式与原书例 7.26 的解析梯度完全一致。对 Black–Scholes 公式，一次前向计算（28 个节点的计算图）加一次反向扫描，就同时得到了 Delta、Vega、Rho 和对到期时间的导数，与解析公式一致到全部有效数字。若用有限差分，每个参数至少要多算一次（中心差分要两次）价格。对一个依赖上百个市场参数（收益率曲线节点、波动率曲面格点）的复杂产品或整个衍生品账簿，AAD 以约 4–5 倍单次定价的代价得到全部敏感度，与参数个数无关——这正是 XVA 和大型账簿风险计算采用 AAD 的原因；长路径蒙特卡洛的计算图过大时，就用检查点技术。分支陷阱也得到了演示：\(f(x)=x-1\) 在 \(x=1\) 处被 AD 算成导数 0。
• D：稀疏 Jacobian。3 次额外求值得到的三对角 Jacobian 与逐列 6 次求值的结果完全相同。\(n\) 再大也只需 3 次。量化中的对应场景是收益率曲线自举/多曲线校准：每个校准工具（存款、FRA、互换）只依赖少数几个曲线节点，Jacobian 呈带状稀疏，着色可以把"每个节点重新定价一遍"的工作量降低一个数量级。

补充：风险贡献就是梯度。组合波动率 \(\sigma_p(w)=\sqrt{w^T\Sigma w}\) 对权重的梯度为 \(\Sigma w/\sigma_p\)，于是 \(w_i\partial\sigma_p/\partial w_i=w_i(\Sigma w)_i/\sigma_p\) 就是第 \(i\) 个资产的风险贡献，且它们之和等于 \(\sigma_p\)（欧拉分解）。对 CVaR、含非线性成本的复杂风险度量，用 AD 框架（PyTorch/JAX）可以自动得到这种分解，不必手推公式。

本章小结

导数可以由有限差分近似、由自动微分精确计算。有限差分要平衡截断误差与舍入误差：前向差分最优步长约 \(\sqrt u\)、精度约 \(\sqrt u\)，中心差分最优步长约 \(u^{1/3}\)、精度约 \(u^{2/3}\)；函数带噪声时要按噪声水平重新选步长并使用公共随机数。Jacobian–向量积、Hessian–向量积只需一次额外的函数或梯度求值；稀疏 Jacobian/Hessian 可用图着色同时估计多列，Hessian 还能利用对称性。自动微分把程序看成计算图：前向模式沿种子方向传播方向导数，适合求 \(Jp\)；反向模式先前向记录、再反向累加伴随变量，梯度代价至多为函数的 4–5 倍且与变量数无关，代价是存储计算图（可用检查点缓解）。前向套反向可以约 12 倍代价得到 Hessian–向量积。AD 的局限在于截断误差的导数和代码分支。

练习

基础

1. （原书 7.1）对中心差分，假设函数值的计算误差不超过 \(uL_f\)、三阶导数有界 \(M\)，写出总误差界，并证明最优步长 \(\epsilon\propto u^{1/3}\)，可达精度 \(\propto u^{2/3}\)。 提示：截断误差 \(\frac{M}{6}\epsilon^2\)，舍入误差 \(\frac{uL_f}{\epsilon}\)。
2. （原书 7.2）推导 Hessian–向量积的中心差分公式，并说明其误差阶。 提示：\(\nabla^2f(x)p\approx[\nabla f(x+\epsilon p)-\nabla f(x-\epsilon p)]/(2\epsilon)\)，误差 \(O(\epsilon^2)\)。
3. （原书 7.3）用 Taylor 展开验证 (7.21)。
4. （原书 7.8）写出加法、指数、正切、幂 \(s^t\)（\(s,t\) 都是变量）的前向模式传播规则。 提示：\(z=s^t\) 时 \(D_pz=ts^{t-1}D_ps+s^t\ln s\,D_pt\)。
5. （原书 7.9）在例 7.26 的计算图上完成整个反向扫描，写出各节点"完成"的顺序。

进阶

6. （原书 7.5）画出 (7.22) 的邻接图，验证"节点 1 一种颜色、其余节点另一种颜色"满足 Coleman–Moré 的着色条件。
7. （原书 7.10）写出乘法与余弦运算的反向模式规则，并与前向模式比较计算一个完整梯度所需的工作量。
8. （原书 7.13）设 \(f(x)=\frac12[x^Tx+(a^Tx)^2]\)。分别统计计算 \(f\)、\(\nabla f\)、\(\nabla^2f\) 以及 \(\nabla^2f(x)p\) 的运算量，说明为什么 Newton–CG 只需 Hessian–向量积是一个巨大优势。
9. 在本章实战 B 中，固定公共随机数，把步长 \(h\) 从 \(10\) 逐步减小到 \(10^{-4}\)，观察 Delta 估计的偏差和标准差如何变化；再对 Gamma（二阶差分）重复这一实验，解释为什么 Gamma 的差分估计对步长更敏感。
10. 扩展实战 C 中的 AD 类，使之支持 numpy 数组的逐元素运算（或改用 jax.grad 若已安装），对一个 50 个资产的组合计算 \(\sigma_p=\sqrt{w^T\Sigma w}\) 的梯度，并验证 \(\sum_iw_i\partial\sigma_p/\partial w_i=\sigma_p\)。

原书推荐习题：7.1（中心差分最优步长，直接用于 Greeks 步长选择）、7.2 与 7.3（Hessian–向量积与纯函数值 Hessian 差分公式）、7.5 与 7.6（稀疏 Hessian 与给定稀疏结构的着色）、7.9 与 7.10（手工执行反向模式，理解 AAD）、7.13（运算量统计）。另有 7.4（Hessian 对角非零时邻接图与关联图的关系）、7.7（跟踪例 7.26 的前向模式；原书题目中的指标与正文不完全一致，做题时以正文为准）、7.11 与 7.12（减、乘、除的二阶前向规则，验证 (7.39)）。

原书对照

页码换算：原书正文页码 = PDF 页码减 20（已用 PDF 页眉核对；第 1–2 章减 21，第 3–11 章减 20，第 12–16 章减 19，第 17 章至附录减 18）。