第 12 章 反向传播的变体

对应原书第 12 章。第 11 章的反向传播是重大突破，但基本算法在实际问题上慢得令人难以接受，原书说可能要几天甚至几周机时。本章先用一个能看清误差曲面的小例子说明为什么慢，再介绍两类加速办法：一类是从观察中得来的启发式改进（动量、可变学习率），另一类是直接借用数值优化的成熟方法（共轭梯度、Levenberg–Marquardt）。

要强调的是：本章所有算法都用反向传播计算导数（从最后一层向第一层算），区别只在于拿到导数后怎样更新权值。为避免混淆，原书把第 11 章的基本算法称为最速下降反向传播（SDBP, steepest descent backpropagation）。Levenberg–Marquardt 对量化读者尤其重要：它是收益率曲线、波动率曲面等模型校准的主力算法。

学习目标

读完本章，你应当能够：

1. 描述多层网络误差曲面的三个特征：曲率变化剧烈、存在多个局部极小、对称性使原点成为鞍点；据此说明初值应取小随机数。
2. 写出动量反向传播（MOBP），从低通滤波的角度解释它，并用特征值分析证明：对二次函数，无论学习率多大，总有动量系数使其稳定。
3. 执行可变学习率反向传播（VLBP）的三条规则，能手算若干步。
4. 实现共轭梯度反向传播所需的线搜索：区间定位 + 黄金分割，并知道何时重置搜索方向。
5. 从平方和函数的 Hessian 推导 Gauss–Newton 与 Levenberg–Marquardt 算法，解释 \(\mu\) 如何在牛顿法与最速下降之间切换。
6. 用 Marquardt 敏感度反传计算网络误差的 Jacobian 矩阵，执行 LMBP 的四个步骤，并知道它的存储瓶颈。
7. 用 LM 校准一条 Nelson–Siegel 收益率曲线，并用 \(\mathbf{J}^T\mathbf{J}\) 估计参数的标准误。

读前导读

这一章在解决什么问题。 第 11 章解决了"梯度怎么算"，本章解决"拿到梯度后怎么走"。基本的反向传播每一步都沿负梯度走固定步长，就像一个只看当下坡度、步幅永远不变的登山者：在平地上挪不动，在窄谷里左右乱撞。本章的改进分两类。一类是经验规则：给步子加惯性（动量），或者走得顺就迈大步、走错就退回来迈小步（可变学习率）。另一类是数值优化的正规军：共轭梯度和 Levenberg–Marquardt（LM）。

LM 对你最有用。如果你做过收益率曲线拟合、期权波动率曲面校准，或者用过 Excel 的规划求解去拟合一个非线性模型，后台很可能就是 LM。它的核心思想你也熟悉：把非线性模型在当前点做一阶近似（像用久期近似债券价格），近似后就是一个普通的线性回归，解正规方程得到一步；再在新的点重新近似、重新回归，如此反复。这叫 Gauss–Newton。LM 再加一个"刹车"参数 \(\mu\)：近似不靠谱时就把步子缩小、方向转向最速下降。加刹车的写法 \((\mathbf{J}^T\mathbf{J}+\mu\mathbf{I})\) 和岭回归的 \((X^TX+\lambda I)\) 是同一个形式。

需要先想起来的数学。

• 二阶泰勒展开与牛顿法。 在 \(x_0\) 附近 \(F(x)\approx F(x_0)+F'(x_0)\Delta x+\tfrac12F''(x_0)\Delta x^2\)，对 \(\Delta x\) 求最小得 \(\Delta x=-F'/F''\)，这就是牛顿法。多元时 \(F''\) 换成 Hessian 矩阵 \(\nabla^2F\)，除法换成乘逆矩阵。金融对应：用久期加凸性近似价格变动。见 第 00 册第 02 章 导数与泰勒展开 和 第 00 册第 05 章 多元微积分与优化。
• Jacobian 矩阵。 一组函数 \(v_1,\dots,v_N\) 对一组参数 \(x_1,\dots,x_n\) 的全部一阶偏导排成 \(N\times n\) 的矩阵，第 \(i\) 行是第 \(i\) 个残差对各参数的敏感度。例：\(v_1=x_1+2x_2\)，\(v_2=x_1x_2\)，在 \((1,1)\) 处 \(\mathbf{J}=\begin{bmatrix}1&2\\1&1\end{bmatrix}\)。在线性回归里 Jacobian 就是（带负号的）设计矩阵 \(X\)。见第 05 章。
• 特征值与线性递推的稳定性。 递推 \(\mathbf{x}_{k+1}=\mathbf{W}\mathbf{x}_k\) 收敛到 0，当且仅当 \(\mathbf{W}\) 所有特征值的模（复数的绝对值）都小于 1，这叫"谱半径小于 1"。一维就是 AR(1) 模型 \(x_{k+1}=\phi x_k\) 平稳要求 \(|\phi|<1\)，你在 CFA 时间序列里见过。见 第 00 册第 06 章 线性代数速成。
• 一元二次方程的根。 \(\lambda^2-b\lambda+c=0\) 的两根之和为 \(b\)、之积为 \(c\)（韦达定理）；判别式 \(b^2-4c<0\) 时两根是共轭复数，模相等，都等于 \(\sqrt c\)。12.2.3 节的动量稳定性证明就靠这一条。
• 矩阵加 \(\mu\mathbf{I}\) 的效果。 \(\mathbf{A}+\mu\mathbf{I}\) 与 \(\mathbf{A}\) 特征向量相同、每个特征值加 \(\mu\)。所以一个接近奇异（最小特征值接近 0）的矩阵加上 \(\mu\mathbf{I}\) 就变得可逆，条件数变小。见第 06 章。

怎么读这一章。 必读 12.1（为什么慢，看懂误差曲面的三个特征）、12.2（动量）、12.5.1–12.5.2（Gauss–Newton 和 LM），以及 12.7.3 的 Nelson–Siegel 校准，那是最直接的实务应用。12.2.3 的稳定性证明第一次可以只看结论（动量系数足够接近 1 总能稳定），有兴趣再回来读讲解框。12.3 可变学习率的手算例和 12.4 共轭梯度的线搜索细节第一次可以跳过。12.5.3 的 Marquardt 敏感度是第 11 章反向传播的直接变形，读懂第 11 章后会很顺。

12.1 反向传播为什么慢

LMS 在学习率不太大时保证收敛到均方误差最小解，因为单层线性网络的均方误差是二次函数：只有一个驻点，Hessian 处处相同，任何方向上的曲率不变，等高线是椭圆。SDBP 是 LMS 的推广（单层线性时两者相同），但多层网络的误差曲面完全不同。

12.1.1 性能曲面例

为了知道最优解在哪里，原书让一个 1-2-1 网络（两层都是 log-sigmoid）去逼近同结构网络在下列参数下的响应：

响应是 \([-2,2]\) 上一条从 0 升到 1 的双台阶曲线。在 \(p=-2,-1.9,\dots,2\)（41 个点）采样，性能指标取 41 点的误差平方和。每次只让两个参数变化、其余固定在最优值，画出误差曲面：

• \(w^1_{1,1}\) 与 \(w^2_{1,1}\)（图 12.3）：曲面显然不是二次的，曲率变化剧烈——平坦区域适合大学习率，高曲率区域需要小学习率，选不出一个对全局都合适的学习率。平坦区在意料之中：sigmoid 输入很大时输出几乎不变。曲面有多个局部极小：全局极小在 \(w^1_{1,1}=10,\ w^2_{1,1}=1\)，位于一条平行于 \(w^1_{1,1}\) 轴的山谷中；另一个局部极小约在 \(w^1_{1,1}=0.88,\ w^2_{1,1}=38.6\)，位于平行于 \(w^2_{1,1}\) 轴的山谷中。
• \(w^1_{1,1}\) 与 \(b^1_1\)（图 12.4）：极小在 \((10,-5)\)；曲面扭曲，有的地方很陡，有的地方极平。从 \(w^1_{1,1}=0,\ b^1_1=-10\) 出发，梯度几乎为零，最速下降实际上停住不动，尽管离极小点很远。
• \(b^1_1\) 与 \(b^1_2\)（图 12.5）：极小在 \((-5,5)\)。这里能看到多层网络的对称性：存在两个误差相同的极小，第二个对应把第一层两个神经元交换（网络"上下翻转"）。正因为这种对称性，原点往往是误差曲面的鞍点。

白话解释：鞍点是梯度为零、但既不是极小也不是极大的点——沿某些方向是谷底，沿另一些方向是山脊，形状像马鞍。为什么对称性让原点成为鞍点？两个隐层神经元互换位置，网络函数不变，所以误差曲面关于"互换"是对称的，两个等价的极小点分处对称面两侧。如果全部权值都是 0，两个神经元完全相同，处于这面"镜子"的正中间。更糟的是，此时两个神经元收到的梯度也完全相同（反向传播里它们的 \(s\) 一样），它们会永远保持相同、无法分化，等于只有一个神经元。小随机初值的作用就是"打破对称"。

初值的启示：不要设为零（原点倾向于是鞍点）；不要设得很大（远离最优点处曲面非常平坦，sigmoid 饱和、梯度消失）。通常取小随机数，既避开原点又不进入平坦区；Nguyen–Widrow 初始化 [NgWi90] 根据 sigmoid 形状和输入范围确定权值大小，并用偏置把各 sigmoid 的中心散布在输入区间内。另外还应尝试多个初值。

12.1.2 收敛例

下面都用批处理（batching）：整个训练集都呈现之后才更新一次，梯度是各样本梯度的平均，比单样本估计准确得多；训练集完备时就是精确梯度。

只调 \(w^1_{1,1}\) 与 \(w^2_{1,1}\) 的两条 SDBP 轨迹（图 12.6）：

• 轨迹 a 最终收敛到最优解，但非常慢：途中先是中等坡度，然后经过一片很平的区域，最后落入坡度极缓的山谷。增大学习率能加快通过平坦区，但进入山谷后会变得不稳定。
• 轨迹 b 被困在另一条山谷中，收敛到局部极小 \(w^1_{1,1}=0.88,\ w^2_{1,1}=38.6\)。

误差随迭代的曲线（图 12.7，对数横轴）是 SDBP 的典型形态：长时间几乎没有进展，然后在很短时间内骤降。图 12.8 用更大的学习率：开始快得多，但进入含极小点的狭窄山谷后开始发散，在山谷两壁之间来回振荡。

由此得到两个改进思路：(1) 可变学习率——平坦处增大、坡度变陡时减小（问题是算法怎么知道自己在哪儿）；(2) 平滑轨迹——对更新量做平均，滤掉振荡。

12.2 启发式改进之一：动量

12.2.1 低通滤波器

先看一阶滤波器

\(w(k)\) 是输入，\(y(k)\) 是输出，\(\gamma\) 称为动量系数（momentum coefficient）。原书用 \(w(k)=1+\sin(2\pi k/16)\) 演示：\(\gamma=0.9\) 时输出的振荡明显小于输入，\(\gamma=0.98\) 时更小；输出的均值始终等于输入的均值，但 \(\gamma\) 越大响应越慢。一句话：减少振荡，同时跟踪平均值。量化读者会认出这就是指数加权移动平均（EWMA）。

12.2.2 动量反向传播（MOBP）

把参数更新量 \(\Delta\mathbf{W}^m(k)=-\alpha\mathbf{s}^m(\mathbf{a}^{m-1})^T\) 当作滤波器的输入：

原书图 12.10：初值和学习率与图 12.8（原本发散）相同，加上 \(\gamma=0.8\) 的动量后算法变得稳定。动量允许使用更大的学习率而保持稳定；当轨迹连续沿同一方向移动时，更新量不断累积，还会加速收敛。名字由来：它让轨迹倾向于沿原方向继续走，\(\gamma\) 越大"惯性"越大。在狭窄山谷中，横向的来回振荡被平均掉，纵向（沿谷底）的一致分量被保留并累积。

推导拆解：把式 12.9 的递推一路展开（记 \(\mathbf{G}(k)=\mathbf{s}^m(\mathbf{a}^{m-1})^T\) 为第 \(k\) 步的梯度）：

\[\Delta\mathbf{W}(k)=-(1-\gamma)\alpha\big[\mathbf{G}(k)+\gamma\mathbf{G}(k-1)+\gamma^2\mathbf{G}(k-2)+\cdots\big]\]

方括号乘 \((1-\gamma)\) 正是梯度的指数加权移动平均，权重 \((1-\gamma)\gamma^j\) 之和为 1。所以动量的每一步 = 学习率 × 近期梯度的 EWMA。在山谷里，横向梯度一正一负交替，加权平均后相互抵消；沿谷底的梯度符号一致，平均后保留。金融上这和用 EWMA 平滑一个噪声很大的交易信号、只在信号持续同向时才大幅调仓是同一思路。

12.2.3 动量总能使最速下降稳定（P12.2）

这是一个漂亮的结果。对二次函数 \(F=\tfrac12\mathbf{x}^T\mathbf{A}\mathbf{x}+\mathbf{d}^T\mathbf{x}+c\)，带动量的最速下降为 \(\Delta\mathbf{x}_k=\gamma\Delta\mathbf{x}_{k-1}-(1-\gamma)\alpha\mathbf{g}_k\)，\(\mathbf{g}_k=\mathbf{A}\mathbf{x}_k+\mathbf{d}\)。代入 \(\Delta\mathbf{x}_k=\mathbf{x}_{k+1}-\mathbf{x}_k\)：

这是二阶差分方程。令 \(\tilde{\mathbf{x}}_k=[\mathbf{x}_{k-1};\mathbf{x}_k]\)，化为一阶系统 \(\tilde{\mathbf{x}}_{k+1}=\mathbf{W}\tilde{\mathbf{x}}_k+\mathbf{v}\)：

稳定当且仅当 \(\mathbf{W}\) 的特征值模都小于 1。设 \(\mathbf{W}\mathbf{z}=\lambda^w\mathbf{z}\)，\(\mathbf{z}=[\mathbf{z}_1;\mathbf{z}_2]\)，第一块给出 \(\mathbf{z}_2=\lambda^w\mathbf{z}_1\)，第二块给出 \(-\gamma\mathbf{z}_1+\mathbf{T}\mathbf{z}_2=\lambda^w\mathbf{z}_2\)。取 \(\mathbf{z}_2\) 为 \(\mathbf{T}\) 的特征向量（特征值 \(\lambda^t\)），得 \([(\lambda^w)^2-\lambda^t\lambda^w+\gamma]\mathbf{z}_2=\mathbf{0}\)，所以 \(\mathbf{T}\) 的每个特征值对应 \(\mathbf{W}\) 的两个特征值

当 \((\lambda^t)^2<4\gamma\) 时它们是共轭复数，模为 \(|\lambda^w|=\sqrt{(\lambda^t)^2/4+(4\gamma-(\lambda^t)^2)/4}=\sqrt\gamma<1\)——稳定。

推导拆解：这段证明分三步，每步一个技巧。

1. 二阶递推化为一阶。 \(\mathbf{x}_{k+1}\) 依赖 \(\mathbf{x}_k\) 和 \(\mathbf{x}_{k-1}\)，像 AR(2)。把相邻两期堆成一个长向量 \(\tilde{\mathbf{x}}_k\)，就变成"状态向量的一阶递推"，矩阵 \(\mathbf{W}\) 的第一块行只是把 \(\mathbf{x}_k\) 原样搬到上半部分。这和把 AR(2) 写成伴随矩阵（companion matrix）形式完全相同。
2. 把大矩阵的特征值化为小方程。 因为 \(\mathbf{T}\) 和 \(\mathbf{A}\) 只差单位矩阵的倍数，二者特征向量相同；沿每个特征向量方向，问题退化为一维的二次方程 \((\lambda^w)^2-\lambda^t\lambda^w+\gamma=0\)。这一步相当于把耦合的多维问题按"主轴"拆成互不干扰的标量问题。
3. 复根的模。 由韦达定理，两根之积等于常数项 \(\gamma\)；共轭复根模相等，所以每个根的模都是 \(\sqrt\gamma\)。这就是为什么只要进入复根区域，收敛速度就完全由 \(\gamma\) 决定、与学习率无关。

白话解释：复根意味着轨迹以螺旋的方式靠近最优点（有振荡），振幅每步乘以 \(\sqrt\gamma\)。\(\gamma=0.81\) 时每步缩为 0.9 倍，约 22 步缩小到十分之一。\(\mathbf{T}\) 与 \(\mathbf{A}\) 特征向量相同，特征值 \(\lambda^t_i=(1+\gamma)-(1-\gamma)\alpha\lambda_i\)，所以只需

在 \(\gamma=1\) 处两边都等于 2。对 \(\gamma\) 求导：右边斜率为 1，左边斜率为 \(1+\alpha\lambda_i>1\)（强极小时 \(\lambda_i>0\)）。所以当 \(\gamma\) 从 1 稍稍减小时，左边下降得比右边快，不等式成立。结论：对二次函数，无论学习率多大，只要动量系数足够接近 1，带动量的最速下降就稳定，且此时特征值的模为 \(\sqrt\gamma\)。 代价是 \(\sqrt\gamma\) 接近 1 时收敛变慢。原书示例：\(F=x_1^2+25x_2^2\)，\(\alpha=0.041\) 时无动量发散（第 09 章），加 \(\gamma=0.2\) 后轨迹稳定。

12.3 启发式改进之二：可变学习率

单层线性网络的误差曲面是二次的，最大稳定学习率是固定的 \(2/\lambda_{\max}\)。多层网络的曲面形状随区域变化，可以在训练中调整学习率，难点是何时调、调多少。原书介绍一种简单的批量方法 [VoMa88]——可变学习率反向传播（VLBP）：

1. 若一次更新后（整个训练集上的）平方误差增加超过设定比例 \(\zeta\)（通常 1%–5%），则丢弃这次更新，学习率乘以 \(\rho\in(0,1)\)，动量系数（若使用）置零；
2. 若平方误差下降，接受更新，学习率乘以 \(\eta>1\)；若 \(\gamma\) 之前被置零则恢复；
3. 若平方误差增加但不超过 \(\zeta\)，接受更新，学习率不变；若 \(\gamma\) 之前被置零则恢复。

原书用 \(\eta=1.05,\ \rho=0.7,\ \zeta=4\%\) 重做前面的例子（图 12.11、12.12）：轨迹沿直线前进、误差持续下降时学习率不断增大；到达狭窄山谷时学习率迅速减小——每当一步会使误差增加超过 4%，就降学习率并去掉动量，使轨迹能急转弯沿山谷走；之后学习率再次增大；接近收敛越过极小点时再次降低。

例（P12.3）——手算三步。 \(F=x_1^2+25x_2^2\)，\(\mathbf{x}_0=[0.5,0.5]^T\)，\(\alpha=0.05,\ \gamma=0.2,\ \eta=1.5,\ \rho=0.5,\ \zeta=5\%\)。\(F(\mathbf{x}_0)=6.5\)，\(\mathbf{g}_0=[1,25]^T\)。

• 第 1 步：\(\Delta\mathbf{x}_0=-0.8(0.05)[1,25]^T=[-0.04,-1]^T\)，试探点 \([0.46,-0.5]^T\)，\(F=6.4616<6.5\)，接受，\(\alpha\to0.075\)。
• 第 2 步：\(\mathbf{g}_1=[0.92,-25]^T\)，\(\Delta\mathbf{x}_1=0.2[-0.04,-1]^T-0.8(0.075)[0.92,-25]^T=[-0.0632,1.3]^T\)，试探点 \([0.3968,0.8]^T\)，\(F=16.157\)，增幅远超 5%，拒绝；\(\alpha\to0.0375\)，\(\gamma\to0\)。
• 第 3 步：\(\Delta\mathbf{x}_2=-0.0375[0.92,-25]^T=[-0.0345,0.9375]^T\)，试探点 \([0.4255,0.4375]^T\)，\(F=4.966<6.4616\)，接受；\(\gamma\) 恢复为 0.2，\(\alpha\to0.05625\)。

相关变体：Jacobs 的 delta-bar-delta [Jaco88] 给每个参数一个独立的学习率，参数变化方向连续几次相同就增大、方向交替就减小；Tollenaere 的 SuperSAB [Toll90] 类似但规则更复杂；Fahlman 的 Quickprop [Fahl88] 假设误差曲面在极小点附近是开口向上的抛物线，且各权值的影响可以分开考虑。今天深度学习常用的 RMSProp、Adam 等"逐参数自适应学习率 + 动量"的优化器，在思想上正是这条路线的延续。

启发式方法的两大缺点：(1) 需要设定好几个参数（\(\zeta,\rho,\eta\) 等），SDBP 只需一个学习率；复杂的方法可能有五六个参数，性能常对它们敏感且依问题而定；(2) 有时在 SDBP 最终能解决的问题上反而不收敛。算法越复杂，这两点越常见。

12.4 数值优化技术之一：共轭梯度反向传播

第 09 章的共轭梯度不需要二阶导数，又有二次终止性。用于多层网络称为 CGBP。算法回顾：\(\mathbf{p}_0=-\mathbf{g}_0\)；\(\mathbf{x}_{k+1}=\mathbf{x}_k+\alpha_k\mathbf{p}_k\)，\(\alpha_k\) 沿方向极小化；\(\mathbf{p}_k=-\mathbf{g}_k+\beta_k\mathbf{p}_{k-1}\)，\(\beta_k\) 取 HS / FR / PR 之一；未收敛则重复。梯度由反向传播（批量）计算。

但它不能照搬，因为性能指标不是二次的：(1) 不能用式 9.31 解析地求步长，需要线搜索；(2) 一般不会在 \(n\) 步内精确收敛，需要重置。

12.4.1 区间定位

先找一个包含极小点的区间。函数比较法 [Scal85]：先算 \(F(\mathbf{x}_0)\)，再算 \(F(\mathbf{x}_0+\varepsilon\mathbf{p}_0)\)，然后依次在 \(2\varepsilon,4\varepsilon,8\varepsilon,\dots\) 处求值（距离每次加倍），直到函数值上升为止。此时极小点被夹在最后三个求值点中首尾两点构成的区间里（倒数第三个点到最后一个点）。仅凭这些求值无法再缩小区间：极小点可能在区间前半段，也可能在后半段。

12.4.2 区间缩减：黄金分割

至少要在区间内求两个内点 \(c<d\) 的值才能缩小区间（一个内点不提供方向信息）：若 \(F(c)>F(d)\)，极小必在 \([c,b]\)；若 \(F(c)<F(d)\)，极小必在 \([a,d]\)（假设区间内只有一个极小）。黄金分割搜索（Golden Section search）巧妙地安排内点，使每次迭代只需一次新的函数求值。设 \(\tau=0.618\)：

\(c_1=a_1+(1-\tau)(b_1-a_1)\)，\(d_1=b_1-(1-\tau)(b_1-a_1)\)。对 \(k=1,2,\dots\)：

• 若 \(F_c<F_d\)：\(a_{k+1}=a_k\)，\(b_{k+1}=d_k\)，\(d_{k+1}=c_k\)，\(c_{k+1}=a_{k+1}+(1-\tau)(b_{k+1}-a_{k+1})\)；\(F_d\leftarrow F_c\)，\(F_c\leftarrow F(c_{k+1})\)；
• 否则：\(a_{k+1}=c_k\)，\(b_{k+1}=b_k\)，\(c_{k+1}=d_k\)，\(d_{k+1}=b_{k+1}-(1-\tau)(b_{k+1}-a_{k+1})\)；\(F_c\leftarrow F_d\)，\(F_d\leftarrow F(d_{k+1})\)；

直到 \(b_{k+1}-a_{k+1}<tol\)。每次区间缩为原来的 0.618 倍。妙处在于旧的内点恰好成为新区间的内点之一，这正是黄金比例 \(\tau^2=1-\tau\) 的性质。

推导拆解：把区间长度归一为 1，内点放在 \(1-\tau\) 和 \(\tau\) 处。若保留左段 \([0,\tau]\)，旧的左内点 \(1-\tau\) 要能直接充当新区间的右内点，即它应在新区间长度的 \(\tau\) 倍处：\(1-\tau=\tau\cdot\tau\)。解 \(\tau^2+\tau-1=0\) 得 \(\tau=(\sqrt5-1)/2\approx0.618\)。这样每轮只有一个新内点需要求值。对神经网络来说，每次"求值"都要对全部样本做一次前向传播，省下一次求值就省下一次全样本计算。

例（P12.4）。 \(F=\tfrac12\mathbf{x}^T\begin{bmatrix}2&1\\1&2\end{bmatrix}\mathbf{x}\)，\(\mathbf{x}_0=[0.8,-0.25]^T\)，沿 \(\mathbf{p}_0=[-1.35,-0.3]^T\) 搜索。区间定位（\(\varepsilon=0.075\)）：\(F(0)=0.5025\)，\(F(0.075)=0.3721\)，\(F(0.15)=0.2678\)，\(F(0.3)=0.1373\)，\(F(0.6)=0.1893\)——上升了，极小在 \([0.15,0.6]\)。黄金分割：\(c_1=0.15+0.382(0.45)=0.3219\)，\(d_1=0.6-0.382(0.45)=0.4281\)，\(F_c=0.1270\)，\(F_d=0.1085\)；\(F_c>F_d\)，新区间 \([0.3219,0.6]\)，\(c_2=0.4281\)，\(d_2=0.4938\)，\(F_d=0.1232\)；此时 \(F_c<F_d\)，新区间 \([0.3219,0.4938]\)，\(d_3=0.4281\)，\(c_3=0.3876\)，\(F_c=0.1094\)……继续下去收敛到第 09 章解析得到的 \(\alpha_0=0.413\)。

12.4.3 重置与表现

二次函数至多 \(n\) 次迭代收敛（\(n\) 为参数个数）；非二次函数一般不会，而理论也没说一个 \(n\) 步周期之后该用什么方向。最简单的做法：每 \(n\) 次迭代把搜索方向重置为负梯度 [Scal85]。

原书图 12.16 显示 CGBP 的迭代次数远少于前述算法。但这有些误导：每次迭代包含区间定位和黄金分割的多次函数求值，单次迭代的计算量大得多。尽管如此，CGBP 已被证明是多层网络最快的批量训练算法之一 [Char92]。

12.5 数值优化技术之二：Levenberg–Marquardt 算法

LM 算法是牛顿法的变体，专为平方和形式的函数设计——这正是均方误差的形式，所以非常适合网络训练，也是一切非线性最小二乘校准问题的标准工具。更系统的理论（信赖域解释、大残差问题）见第 04 册第 10 章。

12.5.1 从牛顿法到 Gauss–Newton

设 \(F\) 是 \(N\) 个函数的平方和：

梯度第 \(j\) 元 \(\partial F/\partial x_j=2\sum_iv_i\,\partial v_i/\partial x_j\)，矩阵形式

\(\mathbf{J}\) 是 Jacobian 矩阵。Hessian 的 \((k,j)\) 元

若 \(\mathbf{S}\) 很小（残差小，或模型接近线性），\(\nabla^2F\approx2\mathbf{J}^T\mathbf{J}\)，代入牛顿法 \(\mathbf{x}_{k+1}=\mathbf{x}_k-\mathbf{A}_k^{-1}\mathbf{g}_k\) 得 Gauss–Newton 法：

优点是不需要二阶导数。对线性最小二乘（\(\mathbf{v}=\mathbf{t}-\mathbf{G}\mathbf{x}\)），\(\mathbf{S}=\mathbf{0}\)，Gauss–Newton 一步就得到正规方程的解（原书 E12.15）。

推导拆解：梯度与 Hessian 的两行公式各用了一个求导法则。

• 梯度（链式法则）：\(v_i^2\) 对 \(x_j\) 求导，外层 \((\cdot)^2\) 给出 \(2v_i\)，内层给出 \(\partial v_i/\partial x_j\)，再对 \(i\) 求和。"对 \(i\) 求和 \(J_{ij}v_i\)"正是 \(\mathbf{J}^T\mathbf{v}\) 的第 \(j\) 个分量。
• Hessian（乘积法则）：再对 \(x_k\) 求导 \(2v_i\,\partial v_i/\partial x_j\)，这是两个因子的乘积，\((uv)'=u'v+uv'\) 给出两项：一项是两个一阶导相乘（组成 \(\mathbf{J}^T\mathbf{J}\)），另一项是残差 \(v_i\) 乘二阶导（组成 \(\mathbf{S}\)）。

白话解释：Gauss–Newton 还有一个更直观的推法，不用 Hessian。把残差在当前点线性化：\(\mathbf{v}(\mathbf{x}+\Delta\mathbf{x})\approx\mathbf{v}+\mathbf{J}\Delta\mathbf{x}\)。最小化 \(\|\mathbf{v}+\mathbf{J}\Delta\mathbf{x}\|^2\) 是一个普通线性回归：把 \(\mathbf{J}\) 当设计矩阵、\(-\mathbf{v}\) 当因变量，正规方程的解就是式 12.28。所以 Gauss–Newton = "线性化 → 跑一次 OLS → 移到新点 → 再线性化"。被丢掉的 \(\mathbf{S}\) 正是线性化忽略的弯曲部分，残差小或模型近似线性时它不重要。

金融直觉：这和用修正久期近似债券价格、再迭代求到期收益率是同一个逻辑。用 Newton 法从价格反解 YTM 时，每一步都用"价格对收益率的一阶导"做线性外推；价格–收益率曲线弯得越厉害（凸性越大）、初始猜测离答案越远，单步近似就越差，需要更多迭代。

12.5.2 Levenberg–Marquardt

Gauss–Newton 的问题是 \(\mathbf{H}=\mathbf{J}^T\mathbf{J}\) 可能奇异或接近奇异。修正为

若 \(\mathbf{H}\) 的特征值、特征向量为 \(\lambda_i,\mathbf{z}_i\)，则 \(\mathbf{G}\mathbf{z}_i=(\lambda_i+\mu)\mathbf{z}_i\)：特征向量不变，特征值都加 \(\mu\)。\(\mu\) 足够大时 \(\mathbf{G}\) 正定可逆。由此得 Levenberg–Marquardt 算法：

关键性质：\(\mu_k\) 很大时，

即学习率为 \(1/(2\mu_k)\) 的最速下降；\(\mu_k\to0\) 时就是 Gauss–Newton。

金融直觉：式 12.32 与岭回归的系数公式 \((X^TX+\lambda I)^{-1}X^Ty\) 形式完全相同：LM 的每一步，就是对线性化后的问题跑一次岭回归。岭惩罚的效果是把系数往 0 收缩，在这里就是"把这一步 \(\Delta\mathbf{x}\) 往 0 收缩"，即迈小步。\(\mu\) 越大，越不相信线性近似、步子越小越保守。好比风险经理对一个只在小幅变动内可靠的 delta 近似设了头寸限额：近似屡次失败（误差不降反升）就收紧限额（增大 \(\mu\)），连续成功就放宽。

\(\mu\) 的调整：初始取小值（如 \(\mu_0=0.01\)）；若一步不能减小 \(F\)，则 \(\mu\) 乘以 \(\vartheta>1\)（如 10）后重做这一步——\(F\) 最终一定会下降，因为这时走的是沿最速下降方向的小步；若一步使 \(F\) 减小，下一步 \(\mu\) 除以 \(\vartheta\)，向 Gauss–Newton 靠拢以加快收敛。这样 LM 兼得牛顿法的速度和最速下降的保证收敛，也就是第 09 章预告的"出现发散时退化为最速下降"的牛顿法变体。

12.5.3 用于多层网络：Jacobian 的计算

若各样本等概率，均方误差正比于误差平方和

误差向量与参数向量为

\(N=Q\times S^M\)，\(n=S^1(R+1)+S^2(S^1+1)+\cdots+S^M(S^{M-1}+1)\)。Jacobian 的每一行对应某个样本的某个输出误差，每一列对应一个权值或偏置。

标准反向传播算的是平方误差对参数的导数 \(\partial(\mathbf{e}_q^T\mathbf{e}_q)/\partial x_l\)；LM 需要的是误差本身的导数 \(\partial e_{k,q}/\partial x_l\)。仿照敏感度，定义 Marquardt 敏感度：

则 Jacobian 元素为

Marquardt 敏感度的反传递推与标准敏感度（式 11.35）完全相同，只有起点不同：

这是一个 \(S^M\times S^M\) 的对角阵——没有标准反传中的因子 2，也没有误差。然后

推导拆解：Marquardt 敏感度与第 11 章敏感度的关系，只差链式法则的最外一层。标准反传求的是 \(\partial(e^2)/\partial n\)；按一元链式法则 \(\dfrac{\partial(e^2)}{\partial n}=2e\cdot\dfrac{\partial e}{\partial n}\)。LM 需要的是去掉外层的 \(\dfrac{\partial e}{\partial n}\)，所以起点去掉了"\(2e\)"这个因子，只剩 \(\partial e/\partial n^M=-\dot f^M(n^M)\)（\(e=t-a\)，对 \(a\) 的导数是 \(-1\)）。之后层与层之间的传递与误差无关，所以递推公式不变。对单输出网络，有 \(\mathbf{s}^m=2e\,\tilde{\mathbf{s}}^m\)：标准敏感度就是 Marquardt 敏感度乘以 \(2e\)。为什么要分开存？因为 LM 需要 Jacobian 的每一行（每个误差各自的导数）来组装 \(\mathbf{J}^T\mathbf{J}\)，而标准反传只给出它们加权求和后的梯度 \(2\mathbf{J}^T\mathbf{v}\)，信息已经被压缩掉了。

注意每个样本要反传 \(S^M\) 个敏感度向量（矩阵的每一列），因为每个样本有 \(S^M\) 个误差，每个误差对应 Jacobian 的一行。

例（P12.5）。 1-1-1 网络，\(f^1(n)=n^2\)，\(f^2(n)=n\)，训练集 \(\{p_1=1,t_1=1\}\)、\(\{p_2=2,t_2=2\}\)，初值 \(W^1=1,b^1=0,W^2=2,b^2=1\)。

• 前向：\(q=1\)：\(n^1=1,a^1=1,a^2=3,e_1=-2\)；\(q=2\)：\(n^1=2,a^1=4,a^2=9,e_2=-7\)。
• Marquardt 敏感度：\(\tilde S^2_1=-1\)，\(\tilde S^1_1=\dot f^1(n^1_1)W^2\tilde S^2_1=2(1)(2)(-1)=-4\)；\(\tilde S^2_2=-1\)，\(\tilde S^1_2=2(2)(2)(-1)=-8\)。
• Jacobian（参数顺序 \(w^1,b^1,w^2,b^2\)）：

12.5.4 LMBP 的四个步骤

1. 把所有输入送入网络，计算输出和误差 \(\mathbf{e}_q=\mathbf{t}_q-\mathbf{a}^M_q\)，按式 12.34 求误差平方和 \(F(\mathbf{x})\)。
2. 计算 Jacobian：用式 12.46 初始化、式 12.47 反传、式 12.48 拼接，再由式 12.43、12.44 得到各元素。
3. 解式 12.32 得 \(\Delta\mathbf{x}_k\)。
4. 用 \(\mathbf{x}_k+\Delta\mathbf{x}_k\) 重算误差平方和。若比第 1 步小，则 \(\mu\leftarrow\mu/\vartheta\)，\(\mathbf{x}_{k+1}=\mathbf{x}_k+\Delta\mathbf{x}_k\)，回到第 1 步；否则 \(\mu\leftarrow\mu\vartheta\)，回到第 3 步。

收敛判据：梯度 \(2\mathbf{J}^T\mathbf{v}\) 的范数小于预定值，或误差平方和降到目标值。

原书图 12.17 画出了第一次迭代所有可能的 LM 步：小 \(\mu\) 时是 Gauss–Newton 方向，大 \(\mu\) 时是最速下降方向，中间值连成一条曲线；\(\mu\) 增大时步子缩短并转向最速下降方向，因此每次迭代总能减小误差。图 12.18（\(\mu_0=0.01\)，\(\vartheta=5\)）的迭代次数比之前所有方法都少。每次迭代的计算量最大（要解 \(n\times n\) 方程），但对中等规模的网络，LMBP 似乎是最快的训练算法 [HaMe94]。

主要缺点：存储。 要存储并求解 \(n\times n\) 的 \(\mathbf{J}^T\mathbf{J}\)，而其他方法只需存 \(n\) 维梯度。参数很多（原书说通常是几千个以上，视内存而定）时不实用。这也是今天的深度网络（参数以百万、十亿计）不用 LM 的原因。

12.6 五种算法的比较

原书此处写"更多变体见第 19 章"，系沿用第一版章号（第二版第 19 章为 ART，与此无关）；训练算法的实际选择见本册第 22 章。

12.7 量化实战

12.7.1 本章在量化中的位置

• LM 是非线性最小二乘校准的主力：收益率曲线拟合（Nelson–Siegel、Svensson）、期权隐含波动率曲面拟合（SVI、SABR）、Heston 等随机波动率模型的参数校准，本质都是"残差平方和"最小化。Jacobian 用解析导数、有限差分或自动微分得到（Marquardt 敏感度就是对网络的解析 Jacobian），\(\mu\) 的自适应保证从粗糙初值出发也能稳健收敛。
• Gauss–Newton 的适用条件：\(\mathbf{S}(\mathbf{x})\) 小，即残差小或模型近似线性时，\(\mathbf{J}^T\mathbf{J}\) 才是好的 Hessian 近似。市场数据噪声大或模型设定不对时残差大，GN 方向可能很差，LM 会自动增大 \(\mu\)。
• \(\mathbf{J}^T\mathbf{J}\) 给出参数的不确定性：在噪声独立同分布的假设下，参数估计的近似协方差为 \(\hat\sigma^2(\mathbf{J}^T\mathbf{J})^{-1}\)，\(\hat\sigma^2=F/(N-n)\)。它可以用来判断校准参数是否可识别——标准误巨大、\(\mathbf{J}^T\mathbf{J}\) 条件数极高，说明数据无法区分某些参数组合，校准结果每天跳来跳去也就不奇怪了。

白话解释：这个公式就是 OLS 系数协方差 \(\hat\sigma^2(X^TX)^{-1}\) 的非线性版本：在最优点附近把模型线性化，Jacobian \(\mathbf{J}\) 扮演设计矩阵 \(X\) 的角色，\(N-n\) 是自由度（样本数减参数数，和 CFA 里回归的 \(n-k-1\) 同一意思）。条件数高相当于回归中的多重共线性：两个参数对残差的影响几乎成比例（Jacobian 的两列几乎平行），数据只能确定它们的某个组合，确定不了各自的值。

• 动量与自适应学习率：训练收益预测网络时用的 SGD+momentum、Adam 等优化器是 MOBP、delta-bar-delta 的后代；"动量 = 低通滤波"的直觉也与交易信号的指数平滑同构。
• 黄金分割线搜索可用于一维参数的高效搜索，但回测目标函数通常噪声大、不是单峰的，盲目优化单一超参数（止损阈值、持仓期）很容易过拟合（第 03 册第 10b 章）。
• 局部极小与对称性：网络训练结果依赖初值，回测中应多次随机初始化取平均或集成，并检验结果的稳定性。

12.7.2 代码一：五种训练算法比较

在原书 12.1 节的问题上（1-2-1、两层 logsig、拟合 12.1–12.2 参数下的网络响应，41 个点）比较 SDBP、MOBP、VLBP、LMBP。全部 7 个参数都参与训练，从 10 组不同的小随机初值出发。梯度用批量反向传播，Jacobian 用 Marquardt 敏感度。

import numpy as np, time
logsig = lambda n: 1 / (1 + np.exp(-n))
S1 = 2
# 参数向量 x 的排列（同 12.36）：[W1(S1), b1(S1), W2(S1), b2]
def unpack(x):
    return x[:S1, None], x[S1:2*S1, None], x[2*S1:3*S1][None, :], x[3*S1]

def forward(x, p):
    W1, b1, W2, b2 = unpack(x)
    a1 = logsig(W1 @ p[None, :] + b1)
    a2 = logsig(W2 @ a1 + b2)                       # 两层都是 logsig（图 12.1）
    return a1, a2.ravel()

def sse_grad(x, p, t):
    """误差平方和及其梯度（批量反向传播）"""
    W1, b1, W2, b2 = unpack(x)
    a1, a2 = forward(x, p); e = t - a2
    s2 = -2 * a2 * (1 - a2) * e                     # (1, Q)
    s1 = a1 * (1 - a1) * (W2.T @ s2[None, :])        # (S1, Q)
    g = np.concatenate([s1 @ p, s1.sum(1), a1 @ s2, [s2.sum()]])
    return e @ e, g

def jacobian(x, p):
    """Marquardt 敏感度 (12.46)-(12.48) 求 J = de/dx，每行对应一个样本"""
    W1, b1, W2, b2 = unpack(x)
    a1, a2 = forward(x, p)
    S2t = -(a2 * (1 - a2))[None, :]                 # 初值 -F'(n^M)，无因子 2 和误差
    S1t = a1 * (1 - a1) * (W2.T @ S2t)
    return np.column_stack([(S1t * p).T, S1t.T, (S2t * a1).T, S2t.T])

# 目标：同结构网络在 (12.1)(12.2) 参数下的响应，41 个点
p = np.linspace(-2, 2, 41)
x_true = np.array([10, 10, -5, 5, 1, 1, -1.0])
t = forward(x_true, p)[1]
x0 = np.zeros(7)

def run(method, iters=5000, alpha=1.0, gamma=0.9):
    x = x0.copy(); F, g = sse_grad(x, p, t); dx = np.zeros_like(x)
    lr, mu = alpha, 0.01
    hist = [F]
    for k in range(iters):
        if method == "SDBP":
            x = x - lr * g
        elif method == "MOBP":
            dx = gamma * dx - (1 - gamma) * lr * g; x = x + dx
        elif method == "VLBP":                      # eta=1.05, rho=0.7, zeta=4%
            dx_try = gamma * dx - (1 - gamma) * lr * g
            F_try, _ = sse_grad(x + dx_try, p, t)
            if F_try > F * 1.04:
                lr *= 0.7; dx = np.zeros_like(x); hist.append(F); continue
            x = x + dx_try; dx = dx_try
            if F_try < F: lr *= 1.05
        elif method == "LMBP":                      # mu0=0.01, theta=10
            e = t - forward(x, p)[1]; J = jacobian(x, p)
            while mu < 1e10:
                step = -np.linalg.solve(J.T @ J + mu * np.eye(7), J.T @ e)
                F_try, _ = sse_grad(x + step, p, t)
                if F_try < F: x = x + step; mu /= 10; break
                mu *= 10
        F, g = sse_grad(x, p, t); hist.append(F)
    return np.array(hist), x

print("方法   学习率  到达SSE<0.05的迭代数(中位)  到达全局极小的初值数/10  10个初值总耗时")
for m, a, iters in [("SDBP", 1, 20000), ("SDBP", 5, 20000), ("MOBP", 5, 20000),
                    ("VLBP", 5, 20000), ("LMBP", None, 300)]:
    ks, n_glob, t0 = [], 0, time.time()
    for seed in range(10):                          # 多初值：原书反复强调
        x0[:] = np.random.default_rng(seed).uniform(-0.5, 0.5, 7)
        h, x = run(m, iters, alpha=a or 1.0)
        if (h < 0.05).any(): ks.append(np.argmax(h < 0.05))
        n_glob += h[-1] < 1e-3
    med = f"{int(np.median(ks))}({len(ks)}个)" if ks else "未到达"
    print(f"{m}  {str(a):>4}  {med:>12}  {n_glob:>14d}  {time.time()-t0:8.1f}s")

# 看 LMBP 从两个初值出发停在哪里：参数顺序 [w1_11, w1_21, b1_1, b1_2, w2_11, w2_12, b2]
for seed in [1, 2]:
    x0[:] = np.random.default_rng(seed).uniform(-0.5, 0.5, 7)
    h, x = run("LMBP", 300)
    print(f"seed {seed}: SSE={h[-1]:.4f}, x={np.round(x, 2)}")

# P12.5：f1=n^2, f2=n 的 1-1-1 网络 Jacobian
W1, b1, W2, b2 = 1.0, 0.0, 2.0, 1.0
J = []
for pq in [1.0, 2.0]:
    n1 = W1 * pq + b1; a1 = n1**2
    S2 = -1.0; S1 = 2 * n1 * W2 * S2
    J.append([S1 * pq, S1, S2 * a1, S2])
print("P12.5 Jacobian:\n", np.array(J))

输出：

方法   学习率  到达SSE<0.05的迭代数(中位)  到达全局极小的初值数/10  10个初值总耗时
SDBP     1       399(8个)               0       2.0s
SDBP     5           未到达               0       2.0s
MOBP     5       43(10个)               1       2.1s
VLBP     5       15(10个)               2       4.1s
LMBP  None        3(10个)               1       0.2s
seed 1: SSE=0.0227, x=[ -0.88  -0.69   0.     0.   -39.05  42.65  -1.8 ]
seed 2: SSE=0.0000, x=[-10.  10.   5.   5.  -1.   1.  -0.]
P12.5 Jacobian:
 [[ -4.  -4.  -1.  -1.]
 [-16.  -8.  -4.  -1.]]

解读。

1. 速度。 走出初始平台（SSE 从约 1.6 降到 0.05 以下）所需的迭代数：SDBP 约 400 步，且有 2 个初值 2 万步都没走出来；MOBP 43 步；VLBP 15 步；LMBP 只要 3 步，总耗时也最少（0.2 秒，对比其他方法 2–4 秒；具体耗时因机器而异）。与原书结论一致。
2. 动量的稳定作用。 学习率 5 时 SDBP 失败（sigmoid 被推入饱和区，梯度消失，停在 SSE 约 4.6 或 11.8 的平台上），而同样学习率的 MOBP 稳定且快——正是 12.2.3 节的理论结论在非二次曲面上的体现。
3. 局部极小是常态。 10 个初值中，最好的算法也只有 1–2 个到达全局极小，其余停在 SSE≈0.023 一类的局部极小。看 seed 1 的终点：第一层权值很小（约 \(\pm0.8\)）、第二层权值很大（约 \(\pm40\)），与原书描述的局部极小（\(w^1_{1,1}=0.88,\ w^2_{1,1}=38.6\)）属于同一类型——用两个几乎线性的隐层神经元的大幅差值去凑台阶。再快的算法也只是更快地掉进最近的坑，多初值是必需的。
4. 对称性。 seed 2 到达的全局极小是 \(\mathbf{W}^1=[-10,10]^T\)、\(\mathbf{b}^1=[5,5]^T\)、\(\mathbf{W}^2=[-1,1]\)、\(b^2=0\)，与真参数 \([10,10,-5,5,1,1,-1]\) 不同，但网络函数完全相同：因为 \(\mathrm{logsig}(-n)=1-\mathrm{logsig}(n)\)，第一个神经元的符号翻转可以被输出权值和偏置补偿。所以比较不同初值的结果应比较误差，而不是比较参数。
5. P12.5 的 Jacobian 与手算一致。

12.7.3 代码二：用 LM 校准 Nelson–Siegel 收益率曲线

Nelson–Siegel 模型 \(y(\tau)=\beta_0+\beta_1\dfrac{1-e^{-\tau/\lambda}}{\tau/\lambda}+\beta_2\Big(\dfrac{1-e^{-\tau/\lambda}}{\tau/\lambda}-e^{-\tau/\lambda}\Big)\) 对 \(\beta\) 线性、对形状参数 \(\lambda\) 非线性，是典型的非线性最小二乘问题。下面用 10 个期限的模拟报价（噪声 3bp）从粗糙初值校准，与 scipy 的 MINPACK 实现对照，并给出参数标准误。

import numpy as np
from scipy.optimize import least_squares
rng = np.random.default_rng(1)

# Nelson-Siegel 收益率曲线：y(tau) = b0 + b1*L1 + b2*L2，tau 为期限(年)，lam 为形状参数
def ns(x, tau):
    b0, b1, b2, lam = x
    u = tau / lam; L1 = (1 - np.exp(-u)) / u; L2 = L1 - np.exp(-u)
    return b0 + b1 * L1 + b2 * L2

tau = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 20, 30.0])
x_true = np.array([4.0, -2.0, 1.5, 1.8])                     # 单位：%
y = ns(x_true, tau) + rng.normal(0, 0.03, tau.size)        # 报价噪声 3bp

def resid(x): return y - ns(x, tau)                          # v(x) = 误差向量
def jac(x, h=1e-7):                                          # 数值 Jacobian dv/dx
    return np.column_stack([(resid(x + h*np.eye(4)[j]) - resid(x - h*np.eye(4)[j])) / (2*h)
                            for j in range(4)])

def levenberg_marquardt(x, mu=0.01, theta=10, tol=1e-8, maxit=200):
    F = resid(x) @ resid(x)
    for k in range(maxit):
        J, v = jac(x), resid(x)
        g = 2 * J.T @ v                                      # (12.22)
        if np.linalg.norm(g) < tol: break                    # 收敛判据：梯度范数
        while mu < 1e10:
            dx = -np.linalg.solve(J.T @ J + mu * np.eye(4), J.T @ v)   # (12.32)
            F_new = resid(x + dx) @ resid(x + dx)
            if F_new < F: x, F, mu = x + dx, F_new, mu / theta; break
            mu *= theta                                      # 失败则向最速下降靠拢
        else:
            break                                            # 已无法再下降（数值极限）
    return x, F, k

x0 = np.array([3.0, -1.0, 0.0, 1.0])                         # 粗糙初值
x_lm, F_lm, k = levenberg_marquardt(x0)
ref = least_squares(resid, x0, method="lm")
print(f"自写 LM: {k} 次迭代, x = {np.round(x_lm, 4)}, SSE = {F_lm:.5f}")
print(f"scipy  : x = {np.round(ref.x, 4)}, SSE = {2*ref.cost:.5f}")
J = jac(x_lm); s2 = F_lm / (tau.size - 4)
se = np.sqrt(np.diag(s2 * np.linalg.inv(J.T @ J)))           # 近似标准误
print("参数近似标准误:", np.round(se, 3))
print("J'J 条件数: %.1e" % np.linalg.cond(J.T @ J))

输出：

自写 LM: 8 次迭代, x = [ 4.0083 -1.9811  1.5154  1.8724], SSE = 0.00316
scipy  : x = [ 4.0083 -1.9811  1.5154  1.8724], SSE = 0.00316
参数近似标准误: [0.028 0.044 0.193 0.166]
J'J 条件数: 1.6e+03

解读。 30 行的 LM 实现 8 次迭代就收敛，结果与 scipy 完全一致。标准误显示：长端水平 \(\beta_0\) 和斜率 \(\beta_1\) 估得很准（标准误 3–4bp），而曲率 \(\beta_2\) 和形状参数 \(\lambda\) 的标准误约 0.2——它们在 10 个期限上的作用相互替代（\(\lambda\) 变大时曲率因子的驼峰右移，可以部分由 \(\beta_2\) 补偿），\(\mathbf{J}^T\mathbf{J}\) 的条件数 1600 反映了这一点。实务中每天重新校准时 \(\beta_2,\lambda\) 跳动大、\(\beta_0,\beta_1\) 稳定，原因就在这里；常见做法是固定 \(\lambda\)（问题变成线性回归）或对其加正则化（第 13a 章）。

本章小结

多层网络的误差曲面不是二次的：曲率变化剧烈，平坦区与狭窄山谷并存，有多个局部极小，对称性使原点成为鞍点，所以初值应取小随机数并多试几组；单一学习率无法同时适应平坦区和山谷，SDBP 的误差曲线呈"长期停滞 + 短时骤降"。启发式改进中，动量是对更新量的一阶低通滤波，平滑振荡、允许更大的学习率，对二次函数可以证明总有动量系数使其稳定；可变学习率按误差变化自适应调整步长并暂停动量，但参数多且敏感。数值优化技术中，CGBP 用区间定位加黄金分割做线搜索、每 \(n\) 步重置；LMBP 用 \(\mathbf{J}^T\mathbf{J}\) 近似 Hessian、用 \(\mu\) 在 Gauss–Newton 与最速下降之间自适应切换，Jacobian 由 Marquardt 敏感度反传得到，对中等规模问题最快，但存储 \(O(n^2)\)。在量化中，LM 是一切非线性最小二乘校准的标准工具，\(\mathbf{J}^T\mathbf{J}\) 还给出参数的可识别性诊断。

练习

基础

1. 单个 logsig 神经元，训练集 \(\{p_1=-3,t_1=0.5\}\)、\(\{p_2=2,t_2=1\}\)，初值 \(w=0.4,\ b=0.15\)（原书 P12.1）。分别求只用第一个样本和用两个样本批处理时的初始负梯度方向。 提示： 第一个样本：\(a=0.2592\)，\(e=0.2408\)，\(s=-0.0925\)，负梯度 \([-sp,-s]=[-0.2774,0.0925]\)；第二个样本：\(a=0.7211\)，\(s=-0.1122\)，负梯度 \([0.2243,0.1122]\)；平均 \([-0.0265,0.1023]\)，与单样本方向差别很大。
2. 证明：对线性网络（\(\mathbf{v}=\mathbf{t}-\mathbf{G}\mathbf{x}\)），\(\mu=0\) 时 LMBP 一步收敛到最小二乘解（原书 E12.15）。 提示： \(\mathbf{J}=-\mathbf{G}\)，\(\Delta\mathbf{x}=(\mathbf{G}^T\mathbf{G})^{-1}\mathbf{G}^T(\mathbf{t}-\mathbf{G}\mathbf{x}_0)\)，故 \(\mathbf{x}_1=(\mathbf{G}^T\mathbf{G})^{-1}\mathbf{G}^T\mathbf{t}\)，与 \(\mathbf{x}_0\) 无关。
3. 继续 P12.4 的黄金分割，再做两次迭代，写出 \(a_k,b_k,c_k,d_k\)。 提示： 当前 \([0.3219,0.4938]\)，\(c=0.3876\ (F=0.1094)\)，\(d=0.4281\ (F=0.1085)\)；\(F_c>F_d\)，取 \([0.3876,0.4938]\)，新 \(d=0.4938-0.382(0.1062)=0.4532\)……区间向 0.413 收缩。
4. 对 \(F=x_1^2+25x_2^2\)、\(\alpha=0.05\)（\(>0.04\)，无动量会发散），用 P12.2 的条件求使带动量最速下降稳定的 \(\gamma\) 范围。 提示： 二次方程 \((\lambda^w)^2-\lambda^t\lambda^w+\gamma=0\) 两根都在单位圆内的充要条件是 \(|\lambda^t|<1+\gamma\)（\(0\le\gamma<1\)），比正文用的复根条件更宽。代入 \(\lambda^t=(1+\gamma)-(1-\gamma)\alpha\lambda\) 得 \(\gamma>\dfrac{\alpha\lambda_{\max}-2}{\alpha\lambda_{\max}+2}\)。这里 \(\alpha\lambda_{\max}=2.5\)，\(\gamma>1/9\approx0.111\)（直接计算 \(\mathbf{W}\) 的谱半径可验证）。

进阶

5. 用 P12.2 的方法分析：对 P9.1 的函数（Hessian 特征值 4 和 16），学习率 \(\alpha=0.2\) 与 \(\alpha=20\) 时分别需要多大的 \(\gamma\) 才稳定（原书 E12.3），并编程画轨迹验证。 提示： 用上题的公式：\(\alpha=0.2\) 时 \(\alpha\lambda_{\max}=3.2\)，\(\gamma>0.231\)；\(\alpha=20\) 时 \(\alpha\lambda_{\max}=320\)，\(\gamma>0.9876\)，此时谱半径接近 1，收敛很慢。
6. 在 12.7.2 节的代码中加入 CGBP（Polak–Ribière + 区间定位 + 黄金分割 + 每 7 步重置），与其他算法比较迭代次数和函数求值次数。 提示： 迭代次数接近 LMBP，但每次迭代的函数求值次数在 10 次以上。
7. 对 1-2-1 网络（logsig 隐层 + 线性输出，初值同第 11 章 11.4.1 节），以 \(p=1\) 和 \(p=0\) 两个样本、目标 \(1+\sin(\pi p/4)\)，手算 LMBP 第一步所需的 Jacobian（原书 E12.14）。 提示： 每个样本一行，7 列；\(\tilde S^2=-1\)，\(\tilde{\mathbf{S}}^1=\mathrm{diag}(a^1(1-a^1))(\mathbf{W}^2)^T(-1)\)；权值列乘以对应输入。
8. 在 12.7.3 节中，把期限点减少为 \(\{1,2,5,10\}\) 年，或把 \(\lambda\) 固定为 1.8 后只估计三个 \(\beta\)，观察标准误和 \(\mathbf{J}^T\mathbf{J}\) 条件数如何变化，并解释。 提示： 期限减少时可识别性变差，\(\beta_2,\lambda\) 的标准误急剧上升；固定 \(\lambda\) 后问题变成线性回归，条件数大幅下降。

原书推荐习题： P12.2、E12.3–E12.6（动量稳定性的特征值分析，理解动量的数学本质）；P12.3、E12.7、E12.8（手算 VLBP）；P12.4、E12.10–E12.13（区间定位与黄金分割）；P12.5、E12.14、E12.15（Jacobian 与 LM 步，线性情形一步收敛）；E12.16（五种算法编程比较，最有实践价值）。

原书对照