量化交易中文教材

第 11 章 非线性方程组

学习目标

读完本章,你应当能够:

  1. 写出求解 \(r(x)=0\) 的牛顿法,证明它在非退化根附近超线性/二次收敛,并解释退化根上只有线性收敛的原因。
  2. 理解非精确牛顿法的强迫参数 \(\eta_k\) 如何控制收敛速度(线性、超线性、二次),以及它与 GMRES、无矩阵 Jacobian–向量积的配合。
  3. 推导 Broyden 更新,说明它为什么是"满足割线方程且改变最小"的秩一修正,并理解初始矩阵 \(B_0\) 在方程组中的重要性。
  4. 用平方和价值函数 \(\tfrac12\|r\|^2\) 做线搜索或信赖域全局化,知道它的伪局部极小只可能出现在 Jacobian 奇异处,并用 Powell 反例说明病态 Jacobian 的危险。
  5. 理解同伦/延拓方法的思想、弧长参数化和预测–校正,以及它在什么情况下仍会失败。
  6. 在量化场景中(风险预算组合、Merton 结构模型、曲线自举)识别非线性方程组,并处理"多个根中只有一个有意义"的问题。

读前导读

这一章在解决什么问题

你最熟悉的非线性方程是求到期收益率(YTM)和内部收益率(IRR):已知债券价格,找一个 \(y\) 使"按 \(y\) 贴现的现金流现值 = 市场价格"。这个方程没有闭式解,计算器和 Excel 的 RATE、IRR、YIELD 都在背后迭代求根,用的正是本章的主角牛顿法:在当前猜测处用切线(也就是久期)近似价格–收益率曲线,看切线在哪里与目标价格相交,作为下一个猜测。

本章把这件事推广到多个方程、多个未知数:\(r(x)=0\),\(r\) 和 \(x\) 都是 \(n\) 维向量。量化里这样的问题很多:曲线自举时让所有报价同时被精确定价、风险平价组合让每个资产的风险贡献等于预算、Merton 模型从股权市值和股权波动率反解资产价值和资产波动率。它们和第 10 章的最小二乘不同:方程必须精确成立(无套利、一致性条件),而不是"尽量接近"。

本章依次回答四个问题:牛顿法为什么在根附近这么快(11.2.1);Jacobian 太贵或系统太大时怎么办(非精确牛顿、Broyden,11.2.2–11.2.3);起点离根很远时怎样防止发散或振荡(价值函数、线搜索、信赖域,11.3);还有更稳的"延拓法"(11.4)。贯穿全章的一个实务教训是:方程可能有多个根,求解器找到的根不一定是你要的那个。非常规现金流的 IRR 可能有两个,就是你见过的例子。

需要先想起来的数学

1. 一元牛顿法。 求 \(g(x)=0\):在 \(x_k\) 处用切线 \(g(x_k)+g'(x_k)(x-x_k)\) 近似 \(g\),令其为零得 \(x_{k+1}=x_k-g(x_k)/g'(x_k)\)。例:求 \(\sqrt2\),即 \(g(x)=x^2-2\),从 \(x_0=1.5\) 出发,\(x_1=1.5-0.25/3=1.41667\),\(x_2=1.414216\),两步就有 6 位有效数字。见 第 00 册第 02 章 导数与泰勒展开。

2. Jacobian 与线性方程组。 多元情形下切线换成"切平面":\(r(x+p)\approx r(x)+J(x)p\),\(J\) 是 \(n\times n\) 的 Jacobian(第 \(i\) 行是第 \(i\) 个方程对各变量的偏导)。牛顿步就是解线性方程组 \(Jp=-r\)。\(J\) 可逆(非奇异,行列式不为零)时有唯一解。见 第 00 册第 05 章 多元微积分与优化 和 第 06 章 线性代数速成。

3. 大 O 与小 o。 \(O(h^2)\) 表示"不超过常数乘以 \(h^2\)";\(o(h)\) 表示"比 \(h\) 小得多",即 \(o(h)/h\to0\)。例:\(h^2=o(h)\),因为 \(h^2/h=h\to0\)。本章用它们描述收敛速度:下一步误差是 \(o(\text{当前误差})\) 叫超线性,是 \(O(\text{当前误差}^2)\) 叫二次。前缀"Q-"只是说这是按相邻两步误差之比(quotient)来定义的。见 第 00 册第 07 章 概率中的分析工具。

4. 范数与条件数。 \(\|v\|\) 是向量长度;矩阵范数 \(\|A\|\) 是"\(A\) 最多把向量拉长多少倍";条件数 \(\kappa(A)=\|A\|\|A^{-1}\|\) 衡量矩阵离"不可逆"有多近,越大越危险。见 第 00 册第 06 章。

5. 积分形式的泰勒定理。 \(r(x+p)-r(x)=\int_0^1J(x+tp)p\,dt\),即"终点值减起点值 = 沿路径把导数累加起来"。和"债券价格变化 = 沿收益率路径对 \(dP/dy\) 积分"是一回事。见 第 00 册第 03 章 积分。

怎么读这一章

核心必读:11.1(特别是"与优化的区别"和"所有根同样好")、11.2.1(牛顿法与定理 11.2)、11.2.3(Broyden)、11.3.1(价值函数的陷阱)、11.5 与 11.6(量化应用与实战)。11.2.2 非精确牛顿只需记住"越接近解,线性方程要解得越准"。11.3.2–11.3.3 的定理证明第一次可跳过,但 Powell 反例的结论(病态 Jacobian 下会悄悄停在错误的点)值得记住。11.2.4 张量方法可以跳过。11.4 延拓法先读 11.4.1 的思想和 11.4.3 的例子,弧长参数化等需要时再看。

建议顺序:11.1 → 11.2.1 → 11.6.1 → 11.2.3 → 11.3.1 → 11.5 → 11.6.2 → 其余。


11.1 问题与特点

很多应用不需要显式优化,而是要找满足一组关系的变量值。当关系是 \(n\) 个等式、\(n\) 个变量时,就是求解非线性方程组

\[ r(x)=0,\qquad r:\mathbb R^n\to\mathbb R^n,\tag{11.1} \]

各分量 \(r_i\) 光滑。满足 (11.1) 的点叫解或根(root)。一般可能无解、唯一解或多解。例如 \(r(x)=(x_2^2-1,\ \sin x_1-x_2)^T\) 有无穷多个解,如 \((3\pi/2,-1)\)、\((\pi/2,1)\)。再如一元方程

\[ r(x)=\sin(5x)-x\tag{11.3} \]

有三个根:\(0\) 和约 \(\pm0.519148\)。

与优化的联系与区别。 牛顿法是两个领域许多算法的核心;线搜索、信赖域、子问题的非精确求解、导数计算、全局收敛等问题在两个领域都重要。有些算法通过最小化 \(\sum r_i^2\) 求解,与第 10 章的最小二乘关系密切,但有两点不同:方程数等于变量数,而且要求所有方程在解处精确成立——方程往往代表守恒律、一致性条件或无套利条件,必须严格满足。此外:

  • 优化要二次收敛需要目标的二阶导数,而方程组只需要一阶导数(Jacobian);
  • 拟牛顿法在方程组中不如在优化中那么有用;
  • 无约束优化中目标函数本身就是自然的价值函数,方程组中则要人为选择价值函数,各有缺点;
  • 与优化中多个局部极小有好坏之分不同,方程组的所有根在数学上同样好。若某个根在物理或经济上没有意义,应当重新表述模型(例如加上正性约束或换变量)。

金融直觉:多 IRR 问题就是"所有根同样好"的现成例子。现金流 \((-100,\ +230,\ -132)\)(先投入、再收回、最后还要付一笔清理费用)的 NPV 方程 \(-100+\dfrac{230}{1+r}-\dfrac{132}{(1+r)^2}=0\),令 \(v=1/(1+r)\) 得 \(132v^2-230v+100=0\),两个根 \(v=\tfrac{10}{11}\) 和 \(v=\tfrac{5}{6}\),即 \(r=10\%\) 和 \(r=20\%\)。两个 IRR 在数学上都对,Excel 返回哪一个取决于你给的初始猜测 guess。CFA 的处理办法是改用 NPV 或 MIRR,也就是原书说的"重新表述模型"。

例 11.1(飞机稳定性):平衡方程 \(F(x)=Ax+\phi(x)=0\),\(A\) 是 \(5\times8\) 常数矩阵,\(\phi\) 是二次项;前 5 个变量是角速度、攻角和侧滑角,后 3 个是控制面偏转。固定控制量就得到 5 个方程 5 个未知数;研究控制变化时的行为,需要对每组控制量解一个系统——一族密切相关的非线性方程组。这正是 11.5 节延拓法的用武之地。

全章假设 \(r\) 在区域 \(D\) 内连续可微,Jacobian \(J(x)=\nabla r(x)\)(第 \(i\) 行为 \(\nabla r_i^T\))存在且连续;部分结果假设 Lipschitz 连续:\(\|J(x_0)-J(x_1)\|\le\beta_L\|x_0-x_1\|\)。\(J(x^*)\) 奇异的根称为退化解(degenerate solution),否则为非退化解。


11.2 局部算法

11.2.1 牛顿法

定理 11.1(Taylor):\(r\) 在凸开集 \(D\) 内连续可微,\(x,x+p\in D\),则

\[ r(x+p)=r(x)+\int_0^1J(x+tp)\,p\,dt.\tag{11.5} \]

于是线性模型 \(M_k(p)=r(x_k)+J(x_k)p\) 的误差是 \(o(\|p\|)\)(\(J\) 连续时)或 \(O(\|p\|^2)\)(\(J\) Lipschitz 时)。

Algorithm 11.1(方程组牛顿法):选 \(x_0\);对 \(k=0,1,\dots\),解牛顿方程

\[ J(x_k)p_k=-r(x_k),\tag{11.9} \]

令 \(x_{k+1}=x_k+p_k\)。

这里用线性模型而不是优化中的二次模型,因为线性模型通常有解且导出快速收敛的算法。几个联系:无约束优化的牛顿法等价于对 \(\nabla f(x)=0\) 用 Algorithm 11.1;原书第 18 章的 SQP 等价于对约束问题的一阶最优性条件(第 12 章的 KKT 条件)用它;\(J(x_k)\) 非奇异时 (11.9) 与 Gauss–Newton 方程 (10.20) 等价。

缺点:起点远离解时行为可能混乱,\(J(x_k)\) 奇异时步甚至无定义;Jacobian 可能难以获得;\(n\) 大时精确求解牛顿方程太贵;根可能退化。退化的例子:\(r(x)=x^2\) 只有退化根 0,牛顿迭代 \(x_{k+1}=x_k-x_k^2/(2x_k)=x_k/2\),即 \(x_k=2^{-k}x_0\),只是线性收敛。

金融直觉:用 YTM 看牛顿法。2 年期、年付息 5% 的债券,面值 100,价格 98。方程是 \(P(y)-98=0\),\(P(y)=\dfrac{5}{1+y}+\dfrac{105}{(1+y)^2}\)。导数 \(P'(y)=-D_{\text{mod}}\cdot P\),就是"负的修正久期乘价格"。所以牛顿步可以写成

\[y_{k+1}=y_k+\frac{P(y_k)-98}{D_{\text{mod}}(y_k)\,P(y_k)},\]
即"价格差 ÷ 美元久期",这正是交易员心算收益率的办法。从 \(y_0=5\%\)(票面利率)出发,迭代得 \(6.0756\%\)、\(6.092281\%\)、\(6.0922847\%\),与真值的误差依次约为 \(1.7\times10^{-4}\)、\(3.9\times10^{-8}\)、\(10^{-15}\),每步有效数字翻倍,这就是二次收敛。快的原因是价格–收益率曲线在根附近很接近它的切线,偏差只来自凸性,是误差的平方量级。 退化根在 IRR 里也有对应:现金流 \((-1,\ +2,\ -1)\) 的 NPV \(=-\big(1-\tfrac1{1+r}\big)^2\) 在 \(r=0\) 处有二重根(NPV 曲线在这里与横轴相切,导数为零)。从 \(r_0=10\%\) 出发,牛顿迭代为 \(4.5\%,\ 2.15\%,\ 1.05\%,\ 0.52\%,\dots\),误差每步只减半,正是上面 \(r(x)=x^2\) 的线性收敛。

定理 11.2:\(r\) 在凸开集 \(D\) 内连续可微,\(x^*\in D\) 为非退化解。则 \(x_k\) 充分接近 \(x^*\) 时,牛顿迭代满足 \(x_{k+1}-x^*=o(\|x_k-x^*\|)\)(局部 Q-超线性);若 \(J\) 在 \(x^*\) 附近 Lipschitz 连续,则 \(x_{k+1}-x^*=O(\|x_k-x^*\|^2)\)(局部 Q-二次)。

证明:\(r(x_k)=r(x_k)-r(x^*)=J(x_k)(x_k-x^*)+o(\|x_k-x^*\|)\)。\(J(x^*)\) 非奇异,所以在 \(x^*\) 的某个球内 \(\|J(x)^{-1}\|\le\beta^*\)。两边左乘 \(J(x_k)^{-1}\):

\[ -p_k=(x_k-x^*)+o(\|x_k-x^*\|)\ \Longrightarrow\ x_{k+1}-x^*=x_k+p_k-x^*=o(\|x_k-x^*\|). \]

Lipschitz 时余项 \(\int_0^1[J(x_k+t(x^*-x_k))-J(x_k)](x_k-x^*)dt\) 的范数是 \(O(\|x_k-x^*\|^2)\),得二次收敛。证毕。

推导拆解:把证明拆成四步,记 \(e_k=x_k-x^*\)(当前误差)。

  1. 用 \(r(x^*)=0\) 凑差:\(r(x_k)=r(x_k)-r(x^*)\)。由 (11.5)(从 \(x^*\) 积到 \(x_k\) 的泰勒定理),这等于 \(J(x_k)e_k\) 加一个余项。\(J\) 连续时余项是 \(o(\|e_k\|)\)。
  2. 两边乘 \(J(x_k)^{-1}\):牛顿步满足 \(p_k=-J(x_k)^{-1}r(x_k)\),所以 \(-p_k=e_k+J(x_k)^{-1}\cdot o(\|e_k\|)\)。\(\|J^{-1}\|\le\beta^*\) 是有界常数,乘上去仍是 \(o(\|e_k\|)\)。这一步是"非退化"假设起作用的地方:若 \(J(x^*)\) 奇异,\(J^{-1}\) 会在根附近爆炸,余项就不再小。
  3. 新误差:\(e_{k+1}=x_k+p_k-x^*=e_k+p_k=o(\|e_k\|)\),即超线性。
  4. 二次的来源:余项是"真实 Jacobian 沿路径的平均"与"\(J(x_k)\)"之差乘以 \(e_k\)。Lipschitz 条件说这个差不超过 \(\beta_L\|e_k\|\),再乘 \(\|e_k\|\) 就是 \(O(\|e_k\|^2)\)。用债券的话说:线性化误差来自凸性,凸性项与 \((\Delta y)^2\) 成正比。

11.2.2 非精确牛顿法

\(n\) 很大时,精确解牛顿方程太贵。非精确牛顿法只要求方向满足

\[ \|r_k+J_kp_k\|\le\eta_k\|r_k\|,\qquad \eta_k\in[0,\eta],\tag{11.18} \]

\(\eta\in[0,1)\) 为常数,\(\eta_k\) 叫强迫参数(forcing parameter;第 06 章 6.2 节称为"强制序列",是同一概念,条件 (11.18) 与第 06 章的 \(\|\nabla^2f_kp_k+\nabla f_k\|\le\eta_k\|\nabla f_k\|\) 形式相同)。收敛理论只依赖 (11.18),最重要的实例是用 GMRES 等 Krylov 子空间方法迭代求解 \(Jp=-r\)(CG 不适用,因为 \(J\) 一般不对称正定)。每次 Krylov 迭代只需一次 Jacobian–向量积 \(Jd\),它不需要显式的 \(J\):有限差分 \(Jd\approx[r(x+\epsilon d)-r(x)]/\epsilon\) 只需一次 \(r\) 求值;自动微分前向模式可以精确计算,代价是 \(r\) 求值的小倍数。GMRES 每步多存一个向量,须周期性重启(常每 10 或 20 步)。

定理 11.3:条件同定理 11.2,\(x_k\) 充分接近 \(x^*\) 时:

  1. 若 \(\eta\) 足够小,Q-线性收敛;
  2. 若 \(\eta_k\to0\),Q-超线性收敛;
  3. 若再有 \(J\) Lipschitz 且 \(\eta_k=O(\|r_k\|)\),Q-二次收敛。

证明要点:把 (11.18) 写成 \(J_kp_k=-r_k+v_k\),\(\|v_k\|\le\eta_k\|r_k\|\),于是 \(\|p_k+J_k^{-1}r_k\|\le\beta^*\eta_k\|r_k\|\)。又 \(\|r(x)\|\le4\|J(x^*)\|\|x-x^*\|\)(在足够小的邻域内),合并得

\[ \|x_k+p_k-x^*\|\le\big(4\|J(x^*)\|\beta^*\eta_k+\beta^*\rho(x_k)\big)\|x_k-x^*\|,\tag{11.23} \]

其中 \(\rho(x)\to0\)(\(x\to x^*\))。取 \(\eta\) 足够小使括号不超过 \(\tfrac12\) 得 (1);\(\eta_k\to0\) 时括号趋于零得 (2)。

实践含义:求解精度应随着接近解而提高——远离解时粗略求解即可,接近解时 \(\eta_k\) 要趋于零,才不会浪费外层的快速收敛。

白话解释:(11.18) 的左边 \(\|r_k+J_kp_k\|\) 是"线性方程 \(J_kp=-r_k\) 没解准的那部分",右边要求它不超过当前方程残差的 \(\eta_k\) 倍。\(\eta_k=0.1\) 就是"线性方程只需解到消除 90% 的残差"。 为什么可以这样偷懒:牛顿法本身是用线性模型近似非线性方程,离根远时线性模型本来就不准,把它解到小数点后十位纯属浪费。好比估值时,WACC 只有两位有效数字,现金流预测就没必要精确到分。但到了根附近,线性模型已经非常准,这时线性方程解不准就会成为瓶颈,所以 \(\eta_k\) 要随 \(\|r_k\|\) 一起减小。 GMRES 是一种迭代求解 \(Jp=-r\) 的方法,每次迭代只需要算一次"\(J\) 乘某个向量",不需要把 \(J\) 整个存下来;这与第 07 章 (7.10) 的 Jacobian–向量积正好配合。

11.2.3 Broyden 方法

割线类(拟牛顿)方法不计算 \(J\),而是维护近似 \(B_k\),使其沿刚走过的步模仿真 Jacobian。步为 \(p_k=-B_k^{-1}r(x_k)\);记 \(s_k=x_{k+1}-x_k\),\(y_k=r(x_{k+1})-r(x_k)\)。由 (11.5),

\[ y_k=\int_0^1J(x_k+ts_k)s_kdt\approx J(x_{k+1})s_k+o(\|s_k\|), \]

所以要求割线方程

\[ B_{k+1}s_k=y_k.\tag{11.26} \]

它是 \(n^2\) 个未知数的 \(n\) 个方程,\(n>1\) 时不能唯一确定 \(B_{k+1}\)(\(n=1\) 时就是标量割线法)。最成功的选择是 Broyden 更新

\[ B_{k+1}=B_k+\frac{(y_k-B_ks_k)s_k^T}{s_k^Ts_k}.\tag{11.27} \]

注意它不要求对称——Jacobian 本来就不对称,这与第 08 章的 Hessian 近似不同。

引理 11.4(Dennis–Schnabel):在所有满足 \(Bs_k=y_k\) 的矩阵中,(11.27) 使 \(\|B-B_k\|_2\) 最小。

证明:对任意满足割线方程的 \(B\),有 \(y_k-B_ks_k=(B-B_k)s_k\),所以

\[ \|B_{k+1}-B_k\|=\Big\|\frac{(B-B_k)s_ks_k^T}{s_k^Ts_k}\Big\|\le\|B-B_k\|\,\Big\|\frac{s_ks_k^T}{s_k^Ts_k}\Big\|=\|B-B_k\|, \]

最后一步用了 \(\|ss^T/s^Ts\|_2=1\)。证毕。

推导拆解:Broyden 更新的两个性质。

  1. 满足割线方程:\(B_{k+1}s_k=B_ks_k+(y_k-B_ks_k)\dfrac{s_k^Ts_k}{s_k^Ts_k}=y_k\)。
  2. 只改动 \(s_k\) 方向:对任何与 \(s_k\) 垂直的向量 \(z\)(\(s_k^Tz=0\)),\(B_{k+1}z=B_kz\)。也就是说,这一步只观测到了沿 \(s_k\) 的信息,就只修正沿 \(s_k\) 的部分,其余方向原样保留。引理 11.4 的"改变最小"就是这个意思。
  3. 引理证明中的不等式用了 \(\|AB\|\le\|A\|\|B\|\);而 \(ss^T/(s^Ts)\) 是"投影到 \(s\) 方向"的矩阵,它把 \(s\) 原样保留、把垂直方向压成零,所以最多把向量"拉长"1 倍,范数为 1。 金融直觉:一维时 Broyden 就是割线法:用最近两次的 \((x,r(x))\) 连一条直线代替切线。解 IRR 时如果不想推导 NPV 的导数,就可以这样做:用两个折现率下的 NPV 连线,看它与零的交点。

Algorithm 11.3(Broyden):选 \(x_0\)、\(B_0\);每步解 \(B_kp_k=-r(x_k)\),沿 \(p_k\) 线搜索得 \(\alpha_k\),\(x_{k+1}=x_k+\alpha_kp_k\),按 (11.27) 更新。

数值例:

\[ r(x)=\begin{bmatrix}(x_1+3)(x_2^3-7)+18\\ \sin(x_2e^{x_1}-1)\end{bmatrix},\tag{11.30} \]

非退化根 \(x^*=(0,1)^T\),起点 \((-0.5,1.4)^T\)(初始误差 0.64),\(B_0=J(x_0)\),单位步长。原书表 11.1 给出:Broyden 法 8 步把误差降到 \(0.87\times10^{-15}\),牛顿法只需 4 步(\(0.12\times10^{-15}\));牛顿法误差的指数每步翻倍(二次),Broyden 的相邻误差比 \(0.097,0.0085,0.47,0.17,0.0033,\dots\) 趋于零但有波动(超线性)。\(\|r(x_k)\|\) 以相似的速度趋于零,因为 \(J(x^*)\) 非奇异时 \(r(x_k)\approx J(x^*)(x_k-x^*)\)。11.6 节的代码会复现这张表。

定理 11.5:在定理 11.2 的假设下,存在 \(\epsilon,\delta>0\),若 \(\|x_0-x^*\|\le\delta\) 且 \(\|B_0-J(x^*)\|\le\epsilon\),则 Broyden 迭代良定且 Q-超线性收敛到 \(x^*\)。

第二个条件在实际中难以保证。与无约束优化不同,\(B_0\) 的选择可能对方程组的性能至关重要,一些实现建议取 \(B_0=J(x_0)\) 或其有限差分近似。即使 \(J\) 稀疏,\(B_k\) 一般也是稠密的;\(n\) 大时可以用有限内存版本(用若干长度为 \(n\) 的向量隐式存储 \(B_k\),以 Sherman–Morrison–Woodbury 公式求解,思路同第 09 章)。

11.2.4 张量方法(简介)

针对退化根(特别是 \(J(x^*)\) 秩为 \(n-1\) 或 \(n-2\) 时),Schnabel–Frank 在牛顿线性模型上加一项二阶张量项

\[ \hat M_k(p)=r(x_k)+J(x_k)p+\tfrac12T_kpp,\tag{11.32} \]

\(T_k\) 不用精确二阶导数(存储约 \(n^3/2\),且模型可能无根),而是选为低秩形式 \(T_kuv=\sum_{j=1}^qa_j(s_{jk}^Tu)(s_{jk}^Tv)\),让模型插值前 \(q\) 个迭代点的函数值(\(q<\sqrt n\),存储 \(2nq\))。若模型无根,就取最小化 \(\|\hat M_k(p)\|\) 的步;若该步不是价值函数的下降方向,就退回牛顿方向。数值测试中标准问题迭代次数一般更少,退化问题上明显优于牛顿法。


11.3 全局化:价值函数、线搜索与信赖域

11.3.1 价值函数及其陷阱

单位步长的牛顿法或 Broyden 法只有在离解足够近时才保证收敛。远离解时,迭代可能发散,也可能循环:\(r(x)=-x^5+x^3+4x\) 的根全部非退化(实根只有 \(0\) 与 \(\pm1.6005\) 三个;原书说"五个根"是按复数域计数,另两个是纯虚根),从 \(x_0=1\) 出发,\(r(1)=4\)、\(r'(1)=2\),于是 \(x_1=-1\);对称地 \(x_2=1\)——牛顿法在 \(\pm1\) 之间永远振荡。

增强稳健性需要价值函数(merit function)来判断候选点是否更接近根。最常用的是平方和

\[ f(x)=\tfrac12\|r(x)\|^2=\tfrac12\sum_ir_i^2(x).\tag{11.34} \]

每个根都是 \(f\) 的全局极小点,但反过来不对:\(f\) 的局部极小点不一定是根。例如 (11.3) 的价值函数除了三个根之外还有许多局部极小(如约 \(\pm1.53053\))。这类伪极小点满足

\[ \nabla f(x^*)=J(x^*)^Tr(x^*)=0,\tag{11.35} \]

而 \(r(x^*)\ne0\),所以只可能出现在 \(J(x^*)\) 奇异的地方。由于这些点会吸引算法,方程组的全局收敛结果不如无约束优化那样令人满意。另一个选择是 \(\ell_1\) 价值函数 \(\|r(x)\|_1\)。

白话解释:价值函数是给"离根还有多远"打的分,分数为 0 就是找到了根。问题在于,分数有局部低谷:在某些点,往任何方向走一点,分数都会变高,但分数并不为 0。(11.35) 说明这种点上 \(J^Tr=0\) 而 \(r\ne0\),即 \(r\) 是一个被 \(J^T\) 映射成零的非零向量,这只有在 \(J\) 奇异时才可能。 一元例子最直观:\(\tfrac12r(x)^2\) 的导数是 \(r(x)r'(x)\)。\(r\ne0\) 时导数为零只能是 \(r'(x)=0\),也就是函数曲线走到了一个"波峰或波谷",但还没碰到横轴。比如 NPV 曲线在某个折现率处出现局部极值却没穿过零,从那附近出发的"按 NPV 绝对值下降"的搜索就会停在那里。

11.3.2 线搜索方法

对 \(f=\tfrac12\|r\|^2\) 用线搜索,方向须是下降方向:

\[ \cos\theta_k=\frac{-p_k^T\nabla f(x_k)}{\|p_k\|\,\|\nabla f(x_k)\|}>0.\tag{11.36} \]

定理 11.6:\(J\) 在水平集 \(\{x:f(x)\le f(x_0)\}\) 的邻域 \(D\) 内 Lipschitz 连续,方向满足 (11.36),步长满足 Wolfe 条件,则 Zoutendijk 条件 \(\sum_k\cos^2\theta_k\|J_k^Tr_k\|^2<\infty\) 成立。

(证明:令 \(\beta_R=\max(\sup_D\|r\|,\sup_D\|J\|)\)——精读笔记指出原书此处把 \(\|r(x)\|\) 误印为 \(f(x)\)——则 \(\|\nabla f(y)-\nabla f(z)\|\le(\beta_L\beta_R+\beta_R^2)\|y-z\|\),\(\nabla f\) Lipschitz 且 \(f\ge0\),套用 Zoutendijk 定理。)

因此只要保证 \(\cos\theta_k\ge\delta>0\),就有 \(J_k^Tr_k\to0\);若再有 \(\|J(x)^{-1}\|\) 在 \(D\) 上有界,则 \(r_k\to0\)。

牛顿方向的 \(\cos\theta_k\)。 \(r_k\ne0\) 时 \(p_k^T\nabla f=-p_k^TJ_k^Tr_k=-\|r_k\|^2<0\),且

\[ \cos\theta_k=\frac{\|r_k\|^2}{\|J_k^{-1}r_k\|\,\|J_k^Tr_k\|}\ge\frac1{\|J_k^T\|\,\|J_k^{-1}\|}=\frac1{\kappa(J_k)}.\tag{11.40} \]

条件数一致有界时没问题;条件数很大时下界接近零,性能可能变差。

推导拆解:(11.40) 的各步。

  1. 牛顿方向 \(p_k=-J_k^{-1}r_k\),价值函数梯度 \(\nabla f=J_k^Tr_k\)。内积 \(p_k^T\nabla f=-r_k^TJ_k^{-T}J_k^Tr_k=-r_k^Tr_k=-\|r_k\|^2\)(\(J^{-T}J^T=I\))。负号说明牛顿方向一定是价值函数的下降方向。
  2. 代入 \(\cos\theta\) 的定义,分母是 \(\|p_k\|\cdot\|\nabla f\|=\|J_k^{-1}r_k\|\cdot\|J_k^Tr_k\|\)。
  3. 用 \(\|J^{-1}r\|\le\|J^{-1}\|\|r\|\) 和 \(\|J^Tr\|\le\|J^T\|\|r\|\),分母不超过 \(\|J^{-1}\|\|J^T\|\|r\|^2\),与分子的 \(\|r\|^2\) 约掉,得 \(1/\kappa(J)\)。 含义:牛顿方向"一定下坡",但条件数大时坡度可以极缓(几乎横着走)。Powell 反例就是让这个夹角逼近 90°,线搜索每步只能挪一点点,最后停在不是根的地方。

Powell 反例:

\[ r(x)=\begin{bmatrix}x_1\\ \dfrac{10x_1}{x_1+0.1}+2x_2^2\end{bmatrix},\qquad J(x)=\begin{bmatrix}1&0\\ \dfrac1{(x_1+0.1)^2}&4x_2\end{bmatrix}.\tag{11.41} \]

唯一解 \(x^*=0\),Jacobian 在整条 \(x_2=0\) 轴上奇异。Powell 证明,在靠近(但不在)\(x_1\) 轴的点处,牛顿步趋于平行于 \(x_2\) 轴,几乎与梯度正交,\(\cos\theta_k\) 可以任意接近零。从 \((3,1)^T\) 出发、用精确线搜索,4 步就收敛到 \((1.8016,0)^T\)——它既不是解,也不是价值函数的驻点(沿 \(-x_1\) 方向两个分量都下降)。教训:病态 Jacobian 下,纯牛顿方向 + 线搜索可能悄无声息地停在错误的地方。

修正方向:

\[ p_k=-(J_k^TJ_k+\tau_kI)^{-1}J_k^Tr_k,\tag{11.42} \]

选 \(\tau_k\) 使 \(\cos\theta_k\ge\delta\);\(\tau_k\to\infty\) 时 \(p_k\) 趋于负梯度方向,总能满足。这正是第 10 章的 LM 方向,可以用 \(\begin{bmatrix}J_k\\\tau^{1/2}I\end{bmatrix}\) 的 QR 分解求,不必形成 \(J^TJ\)。

非精确牛顿方向的 \(\cos\theta_k\)。 把 (11.18) 平方展开可得 \(p^T\nabla f\le\frac{\eta^2-1}2\|r\|^2\),又 \(\|p\|\le\|J^{-1}\|(1+\eta)\|r\|\),于是

\[ \cos\theta_k\ge\frac{1-\eta}{2\kappa(J)}. \]

只要条件数有界,非精确性不损害全局收敛。

Algorithm 11.4(线搜索牛顿法):若牛顿步满足 \(\cos\theta_k\ge\delta\) 则取之,否则用 (11.42) 并选 \(\tau_k\);若 \(\alpha=1\) 满足 Wolfe 条件就取单位步,否则线搜索。

定理 11.7:若该算法用牛顿方向产生 \(x_k\to x^*\),\(r(x^*)=0\),\(J(x^*)\) 非奇异,各 \(r_i\) 二阶可微且 Hessian 有界,单位步满足 Wolfe 条件(\(c_1<\tfrac12\))时即接受,则 Q-二次收敛。原因是单位步最终总被接受,方法退化为纯牛顿法。

11.3.3 信赖域方法

最常用的做法是对 \(f=\tfrac12\|r\|^2\) 用一般信赖域算法,模型 Hessian 取 \(J_k^TJ_k\):

\[ \min_p\ m_k(p)=\tfrac12\|r_k+J_kp\|^2\quad\text{s.t.}\quad\|p\|\le\Delta_k,\tag{11.44} \]
\[ \rho_k=\frac{\|r(x_k)\|^2-\|r(x_k+p_k)\|^2}{\|r(x_k)\|^2-\|r(x_k)+J(x_k)p_k\|^2}.\tag{11.45} \]

Algorithm 11.5:给定 \(\bar\Delta\)、\(\Delta_0\in(0,\bar\Delta)\)、\(\eta\in[0,\tfrac14)\)。每步近似求解 (11.44);若 \(\rho_k<\tfrac14\),\(\Delta_{k+1}=\tfrac14\|p_k\|\);若 \(\rho_k>\tfrac34\) 且 \(\|p_k\|=\Delta_k\),\(\Delta_{k+1}=\min(2\Delta_k,\bar\Delta)\);否则不变。\(\rho_k>\eta\) 时接受该步。

狗腿法:Cauchy 点

\[ p_k^C=-\tau_k\frac{\Delta_k}{\|J_k^Tr_k\|}J_k^Tr_k,\qquad \tau_k=\min\Big\{1,\frac{\|J_k^Tr_k\|^3}{\Delta_k\,r_k^TJ_k(J_k^TJ_k)J_k^Tr_k}\Big\},\tag{11.46–11.47} \]

(模型 Hessian 半正定,无需处理不定情形)与牛顿点 \(p_k^J=-J_k^{-1}r_k\) 之间走折线:若 Cauchy 点已到边界就取它;否则沿 \(p^C\to p^J\) 走到边界或到达 \(p^J\)。它每步只需解一个线性方程组,下降量不少于 Cauchy 点的一个比例。scipy.optimize.root(method='hybr') 调用的 MINPACK hybrd/hybrj 就是 Powell 的混合狗腿法(配合 Broyden 型 Jacobian 更新)。

精确子问题解是 \(p_k=-(J_k^TJ_k+\lambda_kI)^{-1}J_k^Tr_k\)——方程组的 LM 方法在这个意义上是最小二乘 LM 的特例(root(method='lm'))。

定理 11.8–11.9:\(\eta=0\) 时 \(\liminf\|J_k^Tr_k\|=0\);\(\eta\in(0,\tfrac14)\) 且 \(J\) Lipschitz 时 \(\lim\|J_k^Tr_k\|=0\)。

定理 11.10(局部二次收敛):Algorithm 11.5 产生的 \(x_k\) 收敛到非退化解 \(x^*\),\(J\) 在 \(x^*\) 附近 Lipschitz,后期子问题精确求解,则二次收敛。

证明思路(值得体会,它说明"为全局收敛所加的措施不干扰快速局部收敛"):(a) 无论信赖域是否起作用,\(\|p_k\|\le\|J_k^{-1}r_k\|\le\beta^*\|r_k\|\to0\)。(b) 估计 \(|1-\rho_k|\) 的分子:\(r(x_k+p_k)=(r_k+J_kp_k)+w\),\(\|w\|\le(\beta_L/2)\|p_k\|^2\),得分子 \(\le\epsilon(x_k)\|p_k\|^2\),\(\epsilon(x_k)\to0\)。(c) 估计分母:取与 \(p_k\) 等长的牛顿方向步作比较,由子问题最优性得分母 \(\ge\|p_k\|\|r_k\|/\beta^*\)。(d) 合并得 \(|1-\rho_k|\le(\beta^*)^2\epsilon(x_k)\to0\),所以最终 \(\rho_k>\tfrac14\),半径不再缩小;而牛顿步的长度趋于零,最终落在信赖域内被采纳,由定理 11.2 得二次收敛。

另一选择是 \(\ell_1\) 价值函数配 \(\ell_\infty\) 信赖域,\(\min_p\|J_kp+r_k\|_1\) s.t. \(\|p\|_\infty\le\Delta\),可以化为线性规划(Duff–Nocedal–Reid)。


11.4 延拓(同伦)方法

11.4.1 思想

除非 \(J(x)\) 在整个关心区域内都非奇异(常常无法保证),牛顿类方法可能收敛到价值函数的伪极小点。延拓法(continuation,也叫同伦法 homotopy)直接瞄准 \(r(x)=0\) 的解:从一个解显而易见的"容易"方程组出发,逐渐把它变形为原方程组,同时追踪解的移动。标准的同伦映射是

\[ H(x,\lambda)=\lambda r(x)+(1-\lambda)(x-a),\tag{11.59} \]

\(a\in\mathbb R^n\) 固定。\(\lambda=0\) 时 \(H=x-a\),解为 \(x=a\);\(\lambda=1\) 时 \(H=r(x)\)。

朴素做法是让 \(\lambda\) 从 0 小步增加到 1,每个 \(\lambda\) 用牛顿法解 \(H(x,\lambda)=0\)(以上一个解为初值)。满足 \(H(x,\lambda)=0\) 的点集叫零路径(zero path)。若每个 \(\lambda\in[0,1]\) 都有唯一解,这就行得通;但零路径常常会"折回":走到某个 \(\lambda_T\) 后,路径要求 \(\lambda\) 减小才能继续,\(\lambda_T\) 叫转向点(turning point)。朴素做法在转向点处丢失路径。实用的延拓法显式跟踪零路径,允许 \(\lambda\) 时增时减,甚至超出 \([0,1]\)。

11.4.2 弧长参数化与预测–校正

弧长参数化:把 \(x\) 和 \(\lambda\) 都看作路径弧长 \(s\) 的函数,\((x(0),\lambda(0))=(a,0)\)。对 \(H(x(s),\lambda(s))=0\) 求导:

\[ \frac{\partial H}{\partial x}\dot x+\frac{\partial H}{\partial\lambda}\dot\lambda=0,\tag{11.60} \]

切向量 \((\dot x,\dot\lambda)\) 位于 \(n\times(n+1)\) 矩阵 \([\partial_xH\ \ \partial_\lambda H]\) 的零空间中;该矩阵满秩时零空间是一维的,再用归一化 \(\|\dot x\|^2+\dot\lambda^2=1\) 确定长度,用"与上一个切向量夹角小于 \(\pi/2\)"确定方向。切向量可通过带列主元的 QR 分解计算(原书 Procedure 11.7)。于是可以调用标准 ODE 初值问题求解器沿路径积分,直到 \(\lambda(s)=1\)。

预测–校正(predictor–corrector):沿切线走一小步得到预测点 \((x^P,\lambda^P)=(x,\lambda)+\epsilon(\dot x,\dot\lambda)\),它一般不在路径上;再用牛顿迭代把它"拉回"路径。校正时固定最近变化最快的那个分量 \(i\)(可能是 \(\lambda\),也可能是 \(x\) 的某个分量):

\[ \begin{bmatrix}\partial H/\partial x&\partial H/\partial\lambda\\ \multicolumn{2}{c}{e_i^T}\end{bmatrix}\begin{bmatrix}\delta x\\\delta\lambda\end{bmatrix}=\begin{bmatrix}-H\\0\end{bmatrix}. \]

接近 \(\lambda\) 的转向点时,固定 \(\lambda\) 会使系统奇异,必须改为固定 \(x\) 的某个分量——这就是越过转向点的技巧。

白话解释:把零路径想成平面上(一维 \(x\) 时)的一条曲线,横轴 \(\lambda\)、纵轴 \(x\)。朴素做法是"\(\lambda\) 每次往右走一格,求对应的 \(x\)",相当于假设这条曲线是一个函数 \(x(\lambda)\)。但曲线可能像字母 S 一样往回拐:在转向点处它的切线是竖直的,再往右就没有解了。弧长参数化不再把 \(\lambda\) 当作自变量,而是"沿着曲线本身走",曲线往哪拐就跟着往哪拐。 (11.60) 是对恒等式 \(H(x(s),\lambda(s))=0\) 两边关于 \(s\) 求导(链式法则):既然沿路径 \(H\) 始终为零,它的变化率也为零。这给出 \(n\) 个方程,未知数是 \(n+1\) 个分量的切向量,所以切向量只确定到"一个方向",再用长度为 1、方向别掉头两个条件定下来。 金融直觉:你校准一个复杂模型找不到好初值时,常会"从简单情形出发逐步加难度":先校准常数波动率,再逐步打开波动率的期限结构和偏斜,每一步用上一步的结果作初值。这就是朴素延拓。转向点对应"某个参数一调过头,校准就突然找不到解"的情况。

定理 11.11(Watson):\(r\) 二阶连续可微,则对几乎所有 \(a\in\mathbb R^n\),存在从 \((a,0)\) 出发、沿途 \([\partial_xH\ \ \partial_\lambda H]\) 满秩的零路径。若该路径在 \(\lambda\in[0,1)\) 上有界,则它有聚点 \((\bar x,1)\) 满足 \(r(\bar x)=0\);若 \(J(\bar x)\) 非奇异,路径弧长有限。换句话说,除非 \(a\) 选得很不走运,路径要么发散到无穷,要么通向原方程组的解。

11.4.3 失败的例子

例 11.2:\(r(x)=x^2-1\),两个非退化根 \(\pm1\)。取 \(a=-2\):

\[ H(x,\lambda)=\lambda(x^2-1)+(1-\lambda)(x+2)=\lambda x^2+(1-\lambda)x+(2-3\lambda).\tag{11.63} \]

固定 \(\lambda\) 解二次方程,判别式为 \((1-\lambda)^2-4\lambda(2-3\lambda)=13\lambda^2-10\lambda+1\),当

\[ \lambda\in\Big(\frac{5-2\sqrt3}{13},\frac{5+2\sqrt3}{13}\Big)\approx(0.118,\,0.651) \]

时为负,没有实根。从 \((-2,0)\) 出发的路径在 \(\lambda\to0.118\) 时变得无界——定理 11.11 中的"发散"情形。换成 \(a=\tfrac12\),就存在通向 \((1,1)\) 的路径(原书习题 11.12)。

结论:即使是这么简单的系统,延拓法也可能失败;但一般而言它比价值函数方法更可靠,代价是显著更多的函数/导数求值和线性代数运算。


11.5 本章方法在量化中的位置

  • 定价与曲线构建中的求根。 隐含波动率反解是一维求根(牛顿法配 Brent 保护)。从债券或互换报价自举零息曲线,若采用"全局求解"(所有节点同时满足所有报价),就是 \(n\) 个方程 \(n\) 个未知数的非线性方程组,常用牛顿法或 Broyden 法;最终的 Jacobian 还可以直接用于计算各报价的风险敏感度。本章关于 \(B_0\) 的选择、价值函数伪极小、病态 Jacobian(Powell 反例)的讨论,正好对应"曲线构建不收敛"的排查思路。
  • 均衡与一致性条件。 风险平价/风险预算组合的权重由"各资产风险贡献等于给定预算"这一非线性方程组决定;Merton 结构模型中由股权市值和股权波动率反解资产价值和资产波动率是一个 \(2\times2\) 方程组;一般均衡和资产定价模型中的不动点问题也属此类。
  • 多根问题。 金融方程组常常有多个数学上合法的根,但只有一个有经济意义(权重为正、波动率为正、违约概率在 \([0,1]\) 内)。原书的建议是"重新表述模型":换变量(如对数变换保证正性)、把问题写成凸优化的最优性条件,或在线搜索中强制可行性。
  • 延拓思想。 逐步改变相关性、约束强度或风险厌恶系数来追踪最优组合(有效前沿就是一条以风险厌恶为参数的"解路径"),以及从 Black–Scholes 参数出发逐步过渡到复杂模型参数的校准,本质上都是延拓。
  • 大规模。 非精确牛顿 + GMRES 适用于 PDE 定价离散化后的非线性系统,如美式期权的罚函数方法、非线性 HJB 方程。

11.6 量化实战

11.6.1 代码一:复现原书表 11.1,以及牛顿法的循环

import numpy as np

def r(x):
    return np.array([(x[0]+3)*(x[1]**3-7) + 18, np.sin(x[1]*np.exp(x[0]) - 1)])
def J(x):
    c = np.cos(x[1]*np.exp(x[0]) - 1)
    return np.array([[x[1]**3 - 7, 3*x[1]**2*(x[0]+3)],
                     [c*x[1]*np.exp(x[0]), c*np.exp(x[0])]])
xstar, x0 = np.array([0.0, 1.0]), np.array([-0.5, 1.4])

# 牛顿法 Algorithm 11.1(单位步长)
x, err_n = x0.copy(), []
for k in range(4):
    x = x + np.linalg.solve(J(x), -r(x)); err_n.append(np.linalg.norm(x - xstar))
# Broyden 法 (11.24)(11.27),B0 = J(x0),单位步长
x, B, err_b = x0.copy(), J(x0), []
for k in range(8):
    s = np.linalg.solve(B, -r(x)); xn = x + s; y = r(xn) - r(x)
    B = B + np.outer(y - B @ s, s)/(s @ s)            # Broyden 更新 (11.27)
    x = xn; err_b.append(np.linalg.norm(x - xstar))
print("迭代  Broyden ||x-x*||   Newton ||x-x*||")
for k in range(8):
    print(f"{k+1:3d}   {err_b[k]:10.2e}       " + (f"{err_n[k]:10.2e}" if k < 4 else ""))
print("Broyden 相邻误差比:", " ".join(f"{err_b[k+1]/err_b[k]:.3f}" for k in range(7)))

# 牛顿法的循环:r(x) = -x^5 + x^3 + 4x,从 x0 = 1 出发
x = 1.0; seq = [x]
for k in range(4):
    x = x - (-x**5 + x**3 + 4*x)/(-5*x**4 + 3*x**2 + 4); seq.append(x)
print("r=-x^5+x^3+4x 的牛顿迭代:", seq)

输出:

迭代  Broyden ||x-x*||   Newton ||x-x*||
  1     6.20e-02         6.20e-02
  2     5.25e-04         2.11e-04
  3     2.46e-04         1.86e-08
  4     4.30e-05         5.55e-17
  5     1.40e-07       
  6     5.69e-10       
  7     1.76e-12       
  8     8.45e-16       
Broyden 相邻误差比: 0.008 0.469 0.175 0.003 0.004 0.003 0.000
r=-x^5+x^3+4x 的牛顿迭代: [1.0, -1.0, 1.0, -1.0, 1.0]

与原书表 11.1 对照:牛顿法第 2、3 步误差 \(2.11\times10^{-4}\)、\(1.86\times10^{-8}\)(原书 \(0.21\times10^{-3}\)、\(0.18\times10^{-7}\)),每步误差指数翻倍;Broyden 法 8 步到 \(8\times10^{-16}\),相邻误差比 0.008、0.47、0.17、0.003……整体趋于零但有波动,与原书表 11.2 的描述一致(原书首个比值 0.097 是从初始误差 0.64 算起的)。最后一行展示了 11.3.1 节的循环:牛顿法在 \(\pm1\) 之间来回跳,永不收敛。

11.6.2 代码二:风险预算组合与 Merton 模型

风险预算组合。(第 06 章 6.6 节用 Hessian-free Newton–CG 极小化同一问题的凸表述;这里从方程组角度求解,第 17 章 17.7.3 节还会从障碍函数角度再看一次。)给定协方差 \(\Sigma\) 和风险预算 \(b_i>0\)(\(\sum b_i=1\)),要求权重 \(w\) 使每个资产的风险贡献占比 \(w_i(\Sigma w)_i/(w^T\Sigma w)\) 等于 \(b_i\)。一个漂亮的变形是:令 \(y\) 满足

\[ F(y)=\Sigma y-b\oslash y=0,\qquad y>0, \]

(\(\oslash\) 为逐元素除法),则 \(w=y/\mathbf 1^Ty\) 就是所求权重——因为 \(y_i(\Sigma y)_i=b_i\),归一化后风险贡献正比于 \(b_i\)。\(F\) 恰好是严格凸函数 \(\tfrac12y^T\Sigma y-\sum_ib_i\ln y_i\) 的梯度,Jacobian \(\Sigma+\operatorname{diag}(b/y^2)\) 对称正定,所以在 \(y>0\) 上根唯一,牛顿法配回溯(保持 \(y>0\))非常可靠。但 \(F(y)=0\) 在 \(y\) 有负分量的区域还有别的根——这正是"多根中只有一个有意义"的典型。

Merton 模型。 股权是公司资产上的看涨期权:\(E=V\,N(d_1)-De^{-rT}N(d_2)\),且由 Itô 引理 \(\sigma_EE=N(d_1)\sigma_VV\)。已知 \(E\)、\(\sigma_E\)、\(D\),解两个未知数 \(V\)、\(\sigma_V\)。

import numpy as np
from scipy.optimize import root
from scipy.stats import norm

# ---------- (1) 风险预算组合:解 F(y) = Σy - b/y = 0 ----------
rng = np.random.default_rng(11)
n = 50
A = np.column_stack([rng.uniform(0.6, 1.4, n)*0.18,                 # 市场因子(年化波动 18%)
                     rng.standard_normal(n)*0.08,                    # 风格因子
                     rng.standard_normal(n)*0.05])                   # 行业因子
Sigma = A @ A.T + np.diag(rng.uniform(0.15, 0.40, n)**2)             # 年化协方差 = 因子 + 特质
b = np.ones(n)/n; b[:10] = 2.0/n; b /= b.sum()                       # 前 10 只给双倍风险预算

F  = lambda y: Sigma @ y - b/y
JF = lambda y: Sigma + np.diag(b/y**2)        # Jacobian 对称正定:这是凸函数 ½y'Σy - Σb ln y 的梯度
def newton_ls(y, tol=1e-12):
    for k in range(50):
        Fy = F(y)
        if np.linalg.norm(Fy) < tol: return y, k
        p = np.linalg.solve(JF(y), -Fy)
        a = 1.0
        while np.any(y + a*p <= 0) or 0.5*np.sum(F(y+a*p)**2) > (1-1e-4*a)*0.5*np.sum(Fy**2):
            a *= 0.5                           # 回溯:保持 y>0 并使价值函数 ½||F||² 充分下降
        y = y + a*p
    return y, k
y0 = 1/np.sqrt(np.diag(Sigma))                # 逆波动率作为初值
y, it = newton_ls(y0.copy())
w = y/y.sum()
rc = w*(Sigma @ w)/(w @ Sigma @ w)            # 风险贡献占比
print(f"牛顿法: {it} 次迭代, ||F||={np.linalg.norm(F(y)):.1e}")
print("风险贡献/预算 的范围: [%.6f, %.6f]" % ((rc/b).min(), (rc/b).max()))
sol = root(F, y0, method='hybr', jac=JF, tol=1e-12)        # MINPACK:Powell 混合狗腿法
print(f"scipy root[hybr]: success={sol.success}, 求值 {sol.nfev}, ||F||={np.linalg.norm(F(sol.x)):.1e}, "
      f"负分量个数 {np.sum(sol.x < 0)}  ← 也是根,但不是我们要的那个")

def broyden(y, B, tol=1e-10, maxit=500):     # Algorithm 11.3:Broyden 更新 + 价值函数回溯
    Fy = F(y)
    for k in range(maxit):
        if np.linalg.norm(Fy) < tol: return y, k
        p = np.linalg.solve(B, -Fy); a = 1.0
        while np.any(y + a*p <= 0) or np.sum(F(y+a*p)**2) > np.sum(Fy**2):
            a *= 0.5
            if a < 1e-10: return y, -k        # 方向不再下降:失败
        s = a*p; Fn = F(y + s); yv = Fn - Fy
        B = B + np.outer(yv - B @ s, s)/(s @ s)   # (11.27)
        y, Fy = y + s, Fn
    return y, maxit
for name, B0 in [("B0 = J(y0)", JF(y0)), ("B0 = c·I", np.mean(np.diag(JF(y0)))*np.eye(n))]:
    yb, kb = broyden(y0.copy(), B0)
    print(f"Broyden, {name:10s}: 迭代 {kb:4d}, ||F||={np.linalg.norm(F(yb)):.1e}")
print("组合年化波动 %.2f%%,最大/最小权重 %.3f / %.3f" % (100*np.sqrt(w @ Sigma @ w), w.max(), w.min()))

# ---------- (2) Merton 模型:由股权市值与股权波动率反解资产价值与资产波动率 ----------
E, sigE, D, rf, T = 30.0, 0.60, 70.0, 0.03, 1.0        # 股权市值、股权波动率、债务面值
def merton(z):
    V, sV = z
    d1 = (np.log(V/D) + (rf + 0.5*sV**2)*T)/(sV*np.sqrt(T)); d2 = d1 - sV*np.sqrt(T)
    return np.array([V*norm.cdf(d1) - D*np.exp(-rf*T)*norm.cdf(d2) - E,     # 股权 = 资产上的看涨期权
                     norm.cdf(d1)*sV*V - sigE*E])                           # 波动率关系 σ_E E = N(d1) σ_V V
sol = root(merton, [E + D, sigE*E/(E + D)], method='hybr')
V, sV = sol.x
d2 = (np.log(V/D) + (rf - 0.5*sV**2)*T)/(sV*np.sqrt(T))
print(f"Merton: V={V:.3f}, σ_V={sV:.4f}, 违约距离 DD={d2:.3f}, 风险中性违约概率 {norm.cdf(-d2):.4%}, 求值 {sol.nfev}")

输出:

牛顿法: 10 次迭代, ||F||=1.6e-16
风险贡献/预算 的范围: [1.000000, 1.000000]
scipy root[hybr]: success=True, 求值 42, ||F||=5.7e-14, 负分量个数 17  ← 也是根,但不是我们要的那个
Broyden, B0 = J(y0): 迭代   76, ||F||=8.7e-11
Broyden, B0 = c·I  : 迭代   79, ||F||=3.0e-11
组合年化波动 17.01%,最大/最小权重 0.042 / 0.011
Merton: V=97.778, σ_V=0.1881, 违约距离 DD=1.843, 风险中性违约概率 3.2701%, 求值 10

解读如下。

  1. 牛顿法 10 步收敛到二次精度,风险贡献与预算之比全部等于 1。
  2. MINPACK 的 hybr 报告"成功",残差也是 \(10^{-14}\),但解有 17 个负分量——它找到了 \(F(y)=0\) 的另一个根,对应含空头的"风险预算组合",在经济上没有意义。求解器没有错:方程组的所有根在数学上同样好。正确的做法是像 newton_ls 那样在线搜索中强制 \(y>0\),或者直接把问题当作凸函数 \(\tfrac12y^T\Sigma y-\sum b_i\ln y_i\) 的极小化问题来解(其定义域本身就要求 \(y>0\))。
  3. Broyden 法需要 76–79 次迭代,远多于牛顿法的 10 次:\(n=50\) 时,秩一更新要很多步才能"学到"整个 Jacobian。本例中 \(B_0=J(y_0)\) 与 \(B_0=cI\) 差别不大,但原书强调在一般方程组中 \(B_0\) 可能至关重要。Jacobian 便宜可得时(本例就是 \(\Sigma+\operatorname{diag}\)),直接用牛顿法。
  4. Merton 模型解得资产价值 97.8、资产波动率 18.8%,风险中性违约距离 1.84,对应风险中性违约概率约 3.3%。初值取 \(V_0=E+D\)、\(\sigma_{V,0}=\sigma_EE/(E+D)\) 是常用的经验初值,离解足够近,局部方法即可。

本章小结

求解 \(r(x)=0\) 的核心是牛顿法:在非退化根附近超线性/二次收敛,退化根上只有线性收敛。非精确牛顿法用强迫参数 \(\eta_k\) 控制线性系统的求解精度,配合 GMRES 和无矩阵 Jacobian–向量积可处理大规模问题。Broyden 法以满足割线方程且变化最小的秩一修正近似 Jacobian,超线性收敛,但 \(B_0\) 的选择比优化中更重要。全局化通常借助价值函数 \(\tfrac12\|r\|^2\):它的伪极小点只出现在 Jacobian 奇异处,牛顿方向的下降性受条件数控制(Powell 反例),可用 LM 型修正;信赖域(狗腿、LM)全局收敛且不破坏局部二次收敛。延拓法沿同伦零路径弧长跟踪,可越过转向点,更可靠但更贵,路径发散时仍会失败。金融方程组常有多个根,只有一个有经济意义,要通过变量变换、可行性保护或改写为凸优化来锁定它。

概念 公式 / 要点
牛顿法 \(J(x_k)p_k=-r(x_k)\);非退化根二次收敛,退化根线性
非精确牛顿 \(|r_k+J_kp_k|\le\eta_k|r_k|\);\(\eta_k\to0\) 超线性,\(\eta_k=O(|r_k|)\) 二次
Broyden 更新 \(B_{k+1}=B_k+\dfrac{(y_k-B_ks_k)s_k^T}{s_k^Ts_k}\);最小变化,不对称
价值函数 \(f=\tfrac12|r|^2\),\(\nabla f=J^Tr\);伪极小只在 \(J\) 奇异处
牛顿方向下降性 \(\cos\theta\ge1/\kappa(J)\);非精确时 \(\ge(1-\eta)/(2\kappa(J))\)
修正方向 \(p=-(J^TJ+\tau I)^{-1}J^Tr\)(LM 型)
狗腿法 Cauchy 点 → 牛顿点的折线;MINPACK hybr
同伦映射 \(H(x,\lambda)=\lambda r(x)+(1-\lambda)(x-a)\);弧长跟踪越过转向点
风险预算 \(\Sigma y=b\oslash y\),\(y>0\),\(w=y/\mathbf 1^Ty\)

练习

基础

  1. 证明 \(\|ss^T/s^Ts\|_2=1\)。(原书习题 11.1)
  2. 对 \(r(x)=x^q\)(\(q>2\) 为整数),证明牛顿法 Q-线性收敛并求收敛比。(原书习题 11.2) 提示:\(x_{k+1}=x_k-x_k/q=(1-1/q)x_k\)。
  3. 验证 \(r(x)=-x^5+x^3+4x\) 从 \(x_0=1\) 出发的牛顿迭代循环,求出它的五个根并验证都是非退化的。(原书习题 11.3) 提示:\(-x(x^4-x^2-4)=0\),\(x=0\) 或 \(x^2=\frac{1\pm\sqrt{17}}2\);实根为 \(0\)、\(\pm1.6005\),另两个 \(x^2=\frac{1-\sqrt{17}}2<0\) 给出纯虚根。各实根处 \(r'(0)=4\)、\(r'(\pm1.6005)\approx-21.1\),都不为零。
  4. 对 \(r(x)=\sin(5x)-x\),平方和价值函数有无穷多个局部极小点。求出它们满足的方程并估计其位置的通式。(原书习题 11.4) 提示:伪极小满足 \(r'(x)=5\cos5x-1=0\) 且 \(r\ne0\)。
  5. 用 SVD 证明 \(\phi(\lambda)=\|(J^TJ+\lambda I)^{-1}J^Tr\|\) 关于 \(\lambda\) 单调递减(除非 \(J^Tr=0\))。(原书习题 11.5)

进阶

  1. 证明定理 11.3(iii):\(J\) Lipschitz 且 \(\eta_k=O(\|r_k\|)\) 时非精确牛顿法 Q-二次收敛。(原书习题 11.7)
  2. 对例 11.2 取 \(a=\tfrac12\),写出 \(H(x,\lambda)\) 的判别式,证明存在从 \((\tfrac12,0)\) 到 \((1,1)\) 的零路径。(原书习题 11.12)
  3. 用预测–校正法(固定 \(\lambda\) 做校正,步长 0.05)数值跟踪练习 7 的零路径。再对 \(a=-2\) 运行,观察在 \(\lambda\approx0.118\) 附近发生什么。
  4. 在 11.6.2 中把风险预算问题改为对 \(u=\ln y\) 求解 \(G(u)=\Sigma e^u-b\oslash e^u=0\)。说明这一变量变换如何自动消除负根,并比较牛顿法迭代次数。
  5. 用 scipy.optimize.root(method='hybr') 求解 Merton 方程组,把初值改为 \(V_0=E\)、\(\sigma_{V,0}=\sigma_E\),观察是否仍收敛;再把杠杆提高(\(D=95\)),讨论初值的重要性,并联系 11.4 节说明可以如何用延拓(逐步增加 \(D\))获得稳健的初值。

原书推荐习题:11.2、11.3(退化根的线性收敛与牛顿循环,体会局部方法的局限)、11.4(价值函数的伪极小)、11.5(LM 步长随 \(\lambda\) 单调,信赖域与 \(\lambda\) 的对应关系)、11.7(非精确牛顿二次收敛证明)。原书第 11 章共有 11.1–11.12 题。


原书对照

本章内容 原书位置 PDF 页码
11.1 问题、例 11.1、与优化的差异 第 11 章引言,式 (11.1)–(11.4) PDF p.296–301
11.2.1 牛顿法、定理 11.1–11.2 §11.1 Newton's Method for Nonlinear Equations PDF p.301–304
11.2.2 非精确牛顿法、定理 11.3 §11.1 Inexact Newton Methods,式 (11.18)–(11.23) PDF p.304–306
11.2.3 Broyden 方法、表 11.1–11.3、定理 11.5 §11.1 Broyden's Method,式 (11.24)–(11.31) PDF p.306–310
11.2.4 张量方法 §11.1 Tensor Methods,式 (11.32)–(11.33) PDF p.310–312
11.3.1–11.3.2 价值函数、线搜索、Powell 反例 §11.2 Merit Functions;Line Search Methods,式 (11.34)–(11.43) PDF p.312–318
11.3.3 信赖域、狗腿法、定理 11.8–11.10 §11.2 Trust-Region Methods,式 (11.44)–(11.58) PDF p.318–324
11.4 延拓/同伦方法、定理 11.11、例 11.2 §11.3 Continuation/Homotopy Methods,式 (11.59)–(11.63) PDF p.324–330
注释与参考、习题 11.1–11.12 Notes and References;Exercises PDF p.330–333

页码换算:本段原书正文页码约等于 PDF 页码减 20。§11.1、§11.2 内部各小节的 PDF 页码为按篇幅估计的大致位置。