量化交易中文教材

元信息:Jorge Nocedal & Stephen J. Wright《Numerical Optimization》,Springer Series in Operations Research。负责 PDF 第 1–219 页。 注意:本 PDF 实际为 第 1 版(1999,ISBN 0-387-98793-2),不是任务描述中的第 2 版。章节编号与第 2 版有差异(如第 1 版第 9 章为“大规模拟牛顿与部分可分优化”,第 2 版拆分重排;第 1 版无第 2 版中的“无导数优化”等章)。页码换算:原书正文页码 ≈ PDF 页码 − 21(例如原书 p.36 = PDF p.57)。PDF 中夹杂前一读者手写的中文生词批注(如“crew:机组人员”),与原书内容无关,已忽略。


前置部分(PDF p.1–21)

  • 封面、丛书页、版权页(PDF p.1–6):Springer 运筹学丛书(编辑 Peter Glynn、Stephen M. Robinson)。作者单位:Nocedal 在 Northwestern University ECE 系,Wright 在 Argonne National Laboratory 数学与计算机科学部。版权 1999。献词给双方父母。
  • 前言(Preface,PDF p.7–10):
    • 本书目标:全面讲述求解连续优化(continuous optimization)问题的最强大、最前沿的技术;先讲每个算法的动机以培养直觉,数学要求尽量少。
    • 由于聚焦连续问题,省略了离散优化与随机优化;但大量应用可表述为连续优化,例如飞行器/机械臂最优轨迹、地震数据反演、在可接受风险下最大化期望收益的投资组合设计、化工过程控制、汽车/飞机零件形状优化。
    • 强调大规模技术:内点法、非精确牛顿法、有限内存方法、部分可分函数、自动微分;对信赖域方法和序列二次规划(SQP)的处理比已有教材更深入;“核心课程”内容包括约束优化理论、牛顿与拟牛顿法、非线性最小二乘与非线性方程、单纯形法、罚函数与障碍法。
    • 读者:工程、运筹、计算机、数学系研究生,可支撑两学期课程;也面向实践者自学。先修:线性代数(含数值线性代数)和微积分,附录汇总了相关背景。多数章节配有简单编程练习。
    • 写作风格:口语化、重视动机。典型章结构:先非严格讨论(配图),再推动并正式陈述算法,然后严格陈述(多数给出证明)主要理论结果,证明可跳过。
    • 未涵盖:网络优化、整数规划、随机规划、非光滑优化、全局优化;非线性规划内点法和互补问题算法也未详细讨论(第 1 版的局限)。
    • 在线资源:NEOS Guide(含投资组合优化等案例)及软件指南。
    • 结语引 Fletcher:优化是“理论与计算、启发式与严谨的迷人结合”。
  • 目录(PDF p.11–20):全书 18 章 + 附录 A。第 1 章引言;第 2–9 章无约束优化(基础、线搜索、信赖域、共轭梯度、实用牛顿法、导数计算、拟牛顿、大规模拟牛顿与部分可分);第 10 章非线性最小二乘;第 11 章非线性方程;第 12 章约束优化理论;第 13–14 章线性规划(单纯形、内点);第 15 章非线性约束优化算法基础;第 16 章二次规划(含“投资组合优化”例子,原书 p.440);第 17 章罚函数、障碍法与增广拉格朗日;第 18 章 SQP;附录 A 背景知识(分析、拓扑、线性代数、矩阵分解、Sherman–Morrison–Woodbury、浮点误差与条件数)。

第 1 章 引言(Introduction,PDF p.21–30)

1.0 概述与建模(PDF p.22–23)

  • 优化无处不在:航空公司排班最小化成本;投资者构建组合以在不承受过高风险的前提下获得高收益;物理系统趋于最小能量;光走最短时间路径。
  • 使用优化需先确定目标(objective,可用单个数衡量的性能指标,如利润、时间、势能),目标依赖于变量/未知量(variables / unknowns),变量常受约束(constraints,如电子密度、贷款利率不能为负)。
  • 识别目标、变量、约束的过程称为建模(modeling),往往是最重要的一步:模型太简单则无实际洞见,太复杂则难以求解。
  • 没有通用优化算法,需针对问题类型选择算法;选择决定求解快慢乃至能否求出解。
  • 求解后要能判断是否成功:最优性条件(optimality conditions)既用于检验当前点是否是解,未满足时也提示如何改进。还可用敏感性分析(sensitivity analysis)考察解对模型和数据变化的敏感程度。

数学表述(Mathematical Formulation,PDF p.23–24)

记号:\(x\) 为变量向量;\(f\) 为目标函数;\(c\) 为约束函数向量。标准形式:

\[ \min_{x\in\mathbb{R}^n} f(x)\quad \text{s.t.}\quad c_i(x)=0,\ i\in\mathcal{E};\qquad c_i(x)\ge 0,\ i\in\mathcal{I}. \tag{1.1} \]
\(\mathcal{E}\)、\(\mathcal{I}\) 分别为等式、不等式约束的指标集。

例(1.2):\(\min (x_1-2)^2+(x_2-1)^2\),s.t. \(x_1^2-x_2\le 0\),\(x_1+x_2\le 2\)。写成标准形式:\(c_1(x)=-x_1^2+x_2\ge0\),\(c_2(x)=-x_1-x_2+2\ge 0\),\(\mathcal{I}=\{1,2\}\),\(\mathcal{E}=\emptyset\)。图 1.1 展示目标函数的等高线(contours,\(f\) 取常值的点集)、可行域(feasible region,满足全部约束的点集)与最优点 \(x^*\)。

要点:常需变换才能写成 (1.1):变量重新编号;最大化 \(f\) 等价于最小化 \(-f\);\(\le\) 约束乘 \(-1\) 变成 \(\ge\)。好的软件会对用户透明地完成这一转换。

例:运输问题(A Transportation Problem,PDF p.25–26)

化工公司有 2 个工厂 \(F_1,F_2\)、12 个零售店 \(R_1,\dots,R_{12}\)。工厂 \(F_i\) 周产能 \(a_i\) 吨,零售店 \(R_j\) 周需求 \(b_j\) 吨,单位运费 \(c_{ij}\)。变量 \(x_{ij}\) 为从 \(F_i\) 运往 \(R_j\) 的吨数:

\[ \min \sum_{ij} c_{ij}x_{ij}\quad\text{s.t.}\quad \sum_{j=1}^{12}x_{ij}\le a_i\ (i=1,2),\quad \sum_{i=1}^{2}x_{ij}\ge b_j\ (j=1,\dots,12),\quad x_{ij}\ge 0. \tag{1.3–1.4} \]
目标与约束都是线性函数,称为线性规划(linear programming, LP)。实际模型还应包含生产和库存成本。

连续 vs. 离散优化(PDF p.25–26)

  • 若工厂生产拖拉机,\(x_{ij}\) 必须为整数,需加约束 \(x_{ij}\in\mathbb{Z}\),成为整数规划(integer programming)。常见误区:先忽略整数约束求实数解再四舍五入,并不能保证接近最优。
  • 离散优化(discrete optimization):解取自有限集合;连续优化:解取自不可数无穷集合(通常是实向量)。连续问题通常更容易,因为函数光滑性使得在某点的目标和约束信息能推断其邻近点的行为;离散问题中“接近”的点函数值可能差别很大,且可行解太多无法穷举。
  • 同时含连续和整数变量的叫混合整数规划(mixed integer programming)。
  • 连续优化算法在离散优化中也重要:分支定界(branch-and-bound)大部分时间在解 LP 松弛(relaxation),通常用单纯形法(第 13 章)。

约束 vs. 无约束优化(PDF p.27)

  • 可按目标与约束的性质(线性、非线性、凸)、变量数(大/小)、光滑性(可微/不可微)分类;最重要的区分是有无约束,本书据此分为两部分。
  • 无约束问题直接来自应用;有时自然约束可安全忽略(不影响最优解);也可作为约束问题的重构——用罚项(penalization terms)替代约束。
  • 约束可以是简单界 \(0\le x_1\le 100\)、一般线性约束 \(\sum_i x_i\le 1\)(注意这正是组合权重的预算约束形式)、或非线性不等式。
  • 目标与约束全线性→LP(管理科学、运筹学大量使用);至少有一个非线性→非线性规划(nonlinear programming),在物理、工程中自然出现,在管理与经济科学中也日益普遍。

全局 vs. 局部优化(PDF p.27)

  • 最快的算法只求局部解(local solution,邻域内目标最小的可行点),未必是全局解(global solution)。全局解难识别更难定位。
  • 重要特例:凸规划(convex programming)中所有局部解都是全局解;LP 属于凸规划。一般非线性问题可能有非全局的局部解。
  • 本书只顺带讨论全局优化;许多全局优化算法通过求解一系列局部优化问题实现。

随机 vs. 确定性优化(PDF p.28)

  • 模型可能依赖建模时未知的量,如运输问题中的需求 \(b_j\),以及许多经济、金融规划模型依赖未来利率走势和经济表现。
  • 建模者可给出若干情景(scenarios)并赋予概率;随机优化(stochastic optimization)利用这些不确定性量化来优化模型的期望表现。
  • 本书只讨论确定性优化(deterministic optimization,模型完全给定);但很多随机优化算法通过求解若干确定性子问题实现。

优化算法(Optimization Algorithms,PDF p.28–29)

优化算法是迭代的:从初始猜测出发生成一列改进估计直至得到解。不同算法的区别在于从一个迭代点移到下一个的策略;会用到 \(f\)、\(c\) 的值以及可能的一阶、二阶导数;有的会累积历史信息,有的只用当前点的局部信息。好算法应具备三个性质:

  • 稳健性(robustness):在其适用问题类上、对所有合理初值都表现良好;
  • 效率(efficiency):计算时间和存储不过多;
  • 精确性(accuracy):能精确找到解,且对数据误差和舍入误差不过分敏感。

这些目标可能冲突(收敛快的方法在大规模问题上可能存储过多;稳健方法可能最慢),收敛速度与存储、稳健性与速度的权衡是数值优化的核心议题。优化理论既刻画最优点又是大多数算法的基础。

凸性(Convexity,PDF p.29–30)

  • 凸集(convex set):\(S\subseteq\mathbb{R}^n\) 中任意两点连线段都在 \(S\) 内,即 \(\forall x,y\in S,\ \alpha\in[0,1]:\ \alpha x+(1-\alpha)y\in S\)。
  • 凸函数(convex function):定义域为凸集,且
    \[f(\alpha x+(1-\alpha)y)\le \alpha f(x)+(1-\alpha)f(y),\quad \forall \alpha\in[0,1].\]
    几何上:函数图像位于连接 \((x,f(x))\) 与 \((y,f(y))\) 的弦的下方。光滑凸函数在 \(n=1,2\) 时呈碗状,等高线围成凸集。图 1.3 示例:\(f(x)=(x_1-6)^2+\tfrac{1}{25}(x_2-4.5)^4\)。\(-f\) 凸则 \(f\) 凹(concave)。
  • 无约束算法通常只保证收敛到驻点(stationary point,可能是极大、极小或拐点);若 \(f\) 凸,则收敛点为全局极小点。
  • 凸规划定义:目标函数凸;等式约束函数 \(c_i,\ i\in\mathcal{E}\) 线性;不等式约束函数 \(c_i,\ i\in\mathcal{I}\) 凹(因为约束写成 \(c_i\ge 0\),凹函数的上水平集是凸集)。凸性使收敛结论更强。

注释与参考(PDF p.30)

优化起源于变分法和 Euler、Lagrange 的工作;1940 年代线性规划的发展推动了现代优化。“数学规划”(mathematical programming)一词诞生于 1940 年代,那时 “programming” 尚未与计算机软件绑定,原意包括问题表述和算法设计分析。建模参考 Dantzig、Ahuja–Magnanti–Orlin、Fourer–Gay–Kernighan(AMPL)、Winston。

本章要点

  1. 优化 = 目标 + 变量 + 约束;标准形式 (1.1),等式指标集 \(\mathcal{E}\)、不等式 \(\mathcal{I}\)(约定 \(c_i\ge 0\))。
  2. 分类维度:连续/离散、有/无约束、线性/非线性、局部/全局、确定/随机、光滑/非光滑。本书只处理连续、确定、光滑、主要求局部解的问题。
  3. 凸集、凸函数、凹函数、凸规划的定义;凸问题局部解即全局解。
  4. 算法评价三要素:稳健性、效率、精确性,三者常需权衡。

与量化交易的关联

  • 组合优化:均值–方差组合正是“在可接受风险下最大化期望收益”,前言与本章都把它列为典型应用。权重预算约束 \(\sum_i w_i=1\)、不卖空 \(w_i\ge0\)、个股上限 \(w_i\le u_i\) 分别对应等式约束、界约束;目标中协方差二次型凸,因此标准均值–方差问题是凸规划(实际是凸 QP,见原书第 16 章)。
  • 离散 vs. 连续的误区:交易中手数/整手约束(A 股 100 股一手)、持仓股票数上限(基数约束)都是整数约束;“先解连续再四舍五入”可能显著偏离最优,在小资金或高价股时尤为明显,本章的警告直接适用。
  • 随机优化:基于情景的组合优化、CVaR 优化、资产负债管理属于随机规划,本书不涉及,但其子问题往往是确定性 LP/QP。
  • 建模的权衡:模型过简(如忽略交易成本)无洞见,过复杂(如非凸冲击成本)难求解——这是策略研究中的常见取舍。
  • 凸性判断是实战第一步:凸问题可放心用局部算法和成熟求解器(如 CVXPY/MOSEK),非凸问题(如风险平价的某些表述、带基数约束问题)要警惕局部解。

推荐习题

本章无习题。


第 2 章 无约束优化基础(Fundamentals of Unconstrained Optimization,PDF p.31–53)

2.0 问题与例 2.1(PDF p.32–34)

无约束优化:

\[\min_x f(x),\quad x\in\mathbb{R}^n,\ f:\mathbb{R}^n\to\mathbb{R}\ \text{光滑}. \tag{2.1}\]
通常我们对 \(f\) 缺乏全局认识,只知道在一列点 \(x_0,x_1,\dots\) 上的函数值和可能的导数;算法可以自行选择这些点,目标是可靠地找到解且不用过多计算或存储。函数信息往往昂贵,因此偏好不浪费求值的算法。

例 2.1(最小二乘数据拟合):在时刻 \(t_1,\dots,t_m\) 测得信号 \(y_1,\dots,y_m\),用模型

\[\phi(t;x)=x_1+x_2 e^{-(x_3-t)^2/x_4}+x_5\cos(x_6 t)\]
拟合。定义残差(residuals)\(r_j(x)=y_j-\phi(t_j;x)\),求
\[\min_{x\in\mathbb{R}^6} f(x)=r_1^2(x)+\cdots+r_m^2(x). \tag{2.3}\]
这是非线性最小二乘(nonlinear least-squares)问题。说明即使变量很少(\(n=6\)),若 \(m\) 很大(如 \(10^5\)),每次求值 \(f\) 也很昂贵。设最优解约为 \(x^*\approx(1.1,0.01,1.2,1.5,2.0,1.5)\),\(f(x^*)=0.34\),目标非零说明模型不能精确复现所有数据。问题:如何验证 \(x^*\) 是极小点?这引出“解”的定义。

2.1 什么是解(What Is a Solution?,PDF p.34–39)

定义:

  • 全局极小点(global minimizer):\(f(x^*)\le f(x)\) 对所有 \(x\)(或建模者关心的定义域)成立。难以找到,因为算法只知道局部信息,永远无法确定未采样区域没有“深坑”。
  • 局部极小点(local minimizer):存在 \(x^*\) 的邻域 \(\mathcal{N}\)(包含 \(x^*\) 的开集)使 \(f(x^*)\le f(x),\ \forall x\in\mathcal{N}\)。也称弱局部极小点(weak local minimizer)。
  • 严格局部极小点(strict / strong local minimizer):\(f(x^*)<f(x),\ \forall x\in\mathcal{N},x\ne x^*\)。
  • 孤立局部极小点(isolated local minimizer):存在邻域使 \(x^*\) 是其中唯一的局部极小点。

例:常函数 \(f(x)=2\) 每点都是弱局部极小点;\(f(x)=(x-2)^4\) 在 \(x=2\) 有严格局部极小点。 反例(严格但不孤立):\(f(x)=x^4\cos(1/x)+2x^4\),\(f(0)=0\),二阶连续可微,\(x^*=0\) 是严格局部极小点,但附近有一列严格局部极小点 \(x_n\to0\)。结论:孤立 ⇒ 严格,反之不成立。

图 2.2:多局部极小的函数,算法容易“困”在局部极小点;分子构象问题的势能函数可能有数百万个局部极小。凸函数是重要特例:每个局部极小都是全局极小。

识别局部极小(Recognizing a Local Minimum,PDF p.36–38)

若 \(f\) 二阶连续可微,可只通过梯度 \(\nabla f(x^*)\) 和 Hessian \(\nabla^2 f(x^*)\) 判断。

定理 2.1(Taylor 定理):设 \(f:\mathbb{R}^n\to\mathbb{R}\) 连续可微,\(p\in\mathbb{R}^n\),则

\[f(x+p)=f(x)+\nabla f(x+tp)^Tp,\quad \text{某 } t\in(0,1). \tag{2.4}\]
若 \(f\) 二阶连续可微,则
\[\nabla f(x+p)=\nabla f(x)+\int_0^1\nabla^2 f(x+tp)p\,dt, \tag{2.5}\]
\[f(x+p)=f(x)+\nabla f(x)^Tp+\tfrac12 p^T\nabla^2 f(x+tp)p,\quad \text{某 } t\in(0,1). \tag{2.6}\]

定理 2.2(一阶必要条件):若 \(x^*\) 是局部极小点且 \(f\) 在 \(x^*\) 的开邻域内连续可微,则 \(\nabla f(x^*)=0\)。 证明思路(反证):若 \(\nabla f(x^*)\neq0\),取 \(p=-\nabla f(x^*)\),则 \(p^T\nabla f(x^*)=-\|\nabla f(x^*)\|^2<0\);由连续性存在 \(T>0\) 使 \(p^T\nabla f(x^*+tp)<0,\ t\in[0,T]\);由 (2.4),对任意 \(\bar t\in(0,T]\),\(f(x^*+\bar t p)=f(x^*)+\bar t p^T\nabla f(x^*+tp)<f(x^*)\),矛盾。

满足 \(\nabla f(x^*)=0\) 的点称为驻点(stationary point)。

正定:\(p^TBp>0,\ \forall p\neq0\);半正定:\(p^TBp\ge0,\ \forall p\)。

定理 2.3(二阶必要条件):若 \(x^*\) 是局部极小点且 \(\nabla^2 f\) 在其开邻域内连续,则 \(\nabla f(x^*)=0\) 且 \(\nabla^2 f(x^*)\) 半正定。 证明:反设存在 \(p\) 使 \(p^T\nabla^2 f(x^*)p<0\),由连续性在 \([0,T]\) 上保持负;由 (2.6),\(f(x^*+\bar tp)=f(x^*)+\bar t p^T\nabla f(x^*)+\tfrac12\bar t^2p^T\nabla^2 f(x^*+tp)p<f(x^*)\),矛盾。

定理 2.4(二阶充分条件):若 \(\nabla^2 f\) 在 \(x^*\) 开邻域内连续,\(\nabla f(x^*)=0\) 且 \(\nabla^2 f(x^*)\) 正定,则 \(x^*\) 是严格局部极小点。 证明:由连续性取半径 \(r\) 使 \(\nabla^2 f\) 在球 \(\mathcal{D}=\{z:\|z-x^*\|<r\}\) 内正定;对任意 \(0<\|p\|<r\),\(f(x^*+p)=f(x^*)+\tfrac12 p^T\nabla^2 f(z)p>f(x^*)\),\(z=x^*+tp\in\mathcal{D}\)。

注意:充分条件比必要条件强(保证严格极小),但不是必要的。反例:\(f(x)=x^4\),\(x^*=0\) 是严格局部极小点,但 Hessian 为 0 不正定。

定理 2.5(凸函数):\(f\) 凸时,任何局部极小点都是全局极小点;若还可微,则任何驻点都是全局极小点。 证明:(1) 反设存在 \(z\) 使 \(f(z)<f(x^*)\),线段 \(x=\lambda z+(1-\lambda)x^*\),\(\lambda\in(0,1]\) 上由凸性 \(f(x)\le\lambda f(z)+(1-\lambda)f(x^*)<f(x^*)\);\(x^*\) 的任何邻域都含该线段的一段,故 \(x^*\) 非局部极小。(2) 同取 \(z\),

\[\nabla f(x^*)^T(z-x^*)=\lim_{\lambda\downarrow0}\frac{f(x^*+\lambda(z-x^*))-f(x^*)}{\lambda}\le\lim_{\lambda\downarrow0}\frac{\lambda f(z)+(1-\lambda)f(x^*)-f(x^*)}{\lambda}=f(z)-f(x^*)<0,\]
故 \(\nabla f(x^*)\neq0\)。

总结:所有无约束算法本质上都在寻找 \(\nabla f(\cdot)=0\) 的点。

非光滑问题(Nonsmooth Problems,PDF p.39)

  • 本书“光滑”通常指二阶导存在且连续。一般不连续函数无法确定极小点;若由少数光滑片段组成,可分别极小化各片。
  • 连续但在某些点不可微时,可借助次梯度(subgradient / generalized gradient)识别解(参见 Hiriart-Urruty & Lemaréchal)。图 2.3:极小点位于“折点”(kink),一阶导有跳跃,在其附近 \(f\) 的行为难以预测,一处的信息无法推断邻近点。
  • 特殊非光滑函数如 \(f(x)=\|r(x)\|_1\)、\(f(x)=\|r(x)\|_\infty\) 有专用算法(Fletcher 第 14 章)。

2.2 算法概览(Overview of Algorithms,PDF p.40–51)

所有无约束算法需用户给初始点 \(x_0\)(了解应用的用户可给出合理估计,否则任意选取)。生成迭代序列 \(\{x_k\}\),在无法继续改进或已足够精确时终止;用 \(x_k\) 处(可能还有历史点)的信息找 \(f\) 值更低的 \(x_{k+1}\)。非单调算法(nonmonotone)不要求每步下降,但要求每 \(m\) 步下降:\(f(x_k)<f(x_{k-m})\)。

两种策略:线搜索与信赖域(PDF p.40–42)

  • 线搜索(line search):先选方向 \(p_k\),再近似求解一维问题
    \[\min_{\alpha>0} f(x_k+\alpha p_k) \tag{2.9}\]
    得步长 \(\alpha\)。精确求解代价大且不必要,实际只试有限个步长直到松散逼近极小。
  • 信赖域(trust region):构造在 \(x_k\) 附近与 \(f\) 相似的模型 \(m_k\),在 \(x_k\) 周围的区域内极小化模型:
    \[\min_p m_k(x_k+p),\quad x_k+p\ \text{位于信赖域内}. \tag{2.10}\]
    若候选步不能充分降低 \(f\),说明信赖域过大,缩小后重解。信赖域通常为球 \(\|p\|_2\le\Delta\),\(\Delta>0\) 为信赖域半径(trust-region radius);也可用椭球或盒形。模型通常为二次:
    \[m_k(x_k+p)=f_k+p^T\nabla f_k+\tfrac12 p^TB_kp, \tag{2.11}\]
    \(B_k\) 为 Hessian 或其近似,\(m_k\) 与 \(f\) 在 \(x_k\) 处一阶一致。
  • 示例:\(f(x)=10(x_2-x_1^2)^2+(1-x_1)^2\) 在 \(x_k=(0,1)\) 处 \(\nabla f_k=(-2,20)^T\),\(\nabla^2 f_k=\mathrm{diag}(-38,20)\)(Hessian 不定)。图 2.4 显示模型等高线(\(m=1\)、\(m=12\))和两个不同半径的信赖域:每次缩小信赖域后,新的候选步不仅更短,方向通常也会改变——这是与线搜索(固定方向)的本质区别。
  • 两者区别在于选“方向”和“距离”的顺序:线搜索先定方向再定步长 \(\alpha_k\);信赖域先定最大距离 \(\Delta_k\),再在此约束下求最佳方向与步长,不满意则缩小 \(\Delta_k\) 重试。

线搜索的搜索方向(Search Directions for Line Search Methods,PDF p.42–47)

最速下降方向(steepest descent)\(-\nabla f_k\):由 Taylor 展开 \(f(x_k+\alpha p)=f(x_k)+\alpha p^T\nabla f_k+\tfrac12\alpha^2p^T\nabla^2 f(x_k+tp)p\),沿 \(p\) 的变化率为 \(p^T\nabla f_k\)。单位方向中下降最快者解

\[\min_p p^T\nabla f_k\quad\text{s.t. } \|p\|=1. \tag{2.12}\]
由 \(p^T\nabla f_k=\|\nabla f_k\|\cos\theta\),在 \(\theta=\pi\) 时最小,即 \(p=-\nabla f_k/\|\nabla f_k\|\),与等高线正交(图 2.5)。优点:只需梯度;缺点:在困难问题上可能极其缓慢。

下降方向(descent direction):与 \(-\nabla f_k\) 夹角严格小于 \(\pi/2\) 的方向。由 \(f(x_k+\epsilon p_k)=f(x_k)+\epsilon p_k^T\nabla f_k+O(\epsilon^2)\),且 \(p_k^T\nabla f_k=\|p_k\|\|\nabla f_k\|\cos\theta_k<0\),故足够小的正步长必使 \(f\) 下降。

牛顿方向(Newton direction),“也许是最重要的方向”:由二阶模型

\[f(x_k+p)\approx f_k+p^T\nabla f_k+\tfrac12 p^T\nabla^2 f_k p\overset{\text{def}}{=}m_k(p) \tag{2.13}\]
在 \(\nabla^2 f_k\) 正定时令导数为零:
\[p_k^N=-(\nabla^2 f_k)^{-1}\nabla f_k. \tag{2.14}\]

  • 可靠性:(2.13) 与精确展开 (2.6) 的差别仅在于用 \(\nabla^2 f(x_k)\) 代替 \(\nabla^2 f(x_k+tp)\);若 Hessian 足够光滑(Lipschitz),误差仅 \(O(\|p\|^3)\),\(\|p\|\) 小时非常精确。
  • 下降性:\(\nabla^2 f_k\) 正定时 \(\nabla f_k^Tp_k^N=-p_k^{N\,T}\nabla^2 f_kp_k^N\le-\sigma_k\|p_k^N\|^2<0\)。
  • 有“自然”步长 1;多数实现尽量用单位步长,仅当下降不充分时才调整。
  • Hessian 非正定时牛顿方向可能不存在或不是下降方向,需修正(第 6 章)。
  • 局部收敛快(通常二次);缺点是需要显式计算 Hessian,有时繁琐、易错、昂贵。

拟牛顿方向(quasi-Newton):用近似 \(B_k\) 代替真 Hessian,每步更新以吸收新信息,仍可达超线性收敛。由 (2.5),

\[\nabla f(x+p)=\nabla f(x)+\nabla^2 f(x)p+\int_0^1[\nabla^2 f(x+tp)-\nabla^2 f(x)]p\,dt,\]
积分项为 \(o(\|p\|)\)。令 \(x=x_k\), \(p=x_{k+1}-x_k\) 得 \(\nabla f_{k+1}=\nabla f_k+\nabla^2 f_{k+1}(x_{k+1}-x_k)+o(\|x_{k+1}-x_k\|)\),在解附近(Hessian 正定区)有
\[\nabla^2 f_{k+1}(x_{k+1}-x_k)\approx\nabla f_{k+1}-\nabla f_k. \tag{2.15}\]
要求 \(B_{k+1}\) 满足割线方程(secant equation):
\[B_{k+1}s_k=y_k,\quad s_k=x_{k+1}-x_k,\ y_k=\nabla f_{k+1}-\nabla f_k. \tag{2.16}\]
额外要求:对称、\(B_{k+1}-B_k\) 低秩;\(B_0\) 由用户选。两个最常用公式:

  • SR1(对称秩一):
    \[B_{k+1}=B_k+\frac{(y_k-B_ks_k)(y_k-B_ks_k)^T}{(y_k-B_ks_k)^Ts_k}. \tag{2.17}\]
  • BFGS(Broyden、Fletcher、Goldfarb、Shanno):
    \[B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}. \tag{2.18}\]
    SR1 为秩一修正,BFGS 为秩二修正;两者都满足割线方程并保持对称。若 \(B_0\) 正定且 \(s_k^Ty_k>0\),BFGS 产生的 \(B_k\) 保持正定。方向 \(p_k=-B_k^{-1}\nabla f_k\)(2.19)。
  • 为避免每步分解 \(B_k\),可直接更新逆 \(H_k=B_k^{-1}\),BFGS 的逆形式:
    \[H_{k+1}=(I-\rho_ks_ky_k^T)H_k(I-\rho_ky_ks_k^T)+\rho_ks_ks_k^T,\quad \rho_k=\frac{1}{y_k^Ts_k}. \tag{2.20}\]
    (原书此处说“(2.17) 和 (2.18) 的等价公式”,实际 (2.20) 是 BFGS 的逆更新。)这样 \(p_k=-H_k\nabla f_k\) 只需矩阵–向量乘法,比分解/回代简单。大规模变体(部分可分、有限内存)见第 9 章。

非线性共轭梯度方向(nonlinear conjugate gradient):

\[p_k=-\nabla f(x_k)+\beta_kp_{k-1},\]
\(\beta_k\) 使 \(p_k\) 与 \(p_{k-1}\) 共轭。CG 原为求解对称正定线性方程组 \(Ax=b\),等价于极小化凸二次函数 \(\phi(x)=\tfrac12x^TAx-b^Tx\)(原书此处印作 \(+b^Tx\),为笔误)。非线性 CG 比最速下降有效得多、计算几乎同样简单,不及牛顿/拟牛顿快,但无需存储矩阵(第 5 章)。

上述方向都可直接用于线搜索框架;除 CG 外都有信赖域对应物。

信赖域的模型(Models for Trust-Region Methods,PDF p.47–48)

  • 取 \(B_k=0\) 且欧氏范数信赖域:\(\min_p f_k+p^T\nabla f_k\) s.t. \(\|p\|_2\le\Delta_k\),闭式解 \(p_k=-\Delta_k\nabla f_k/\|\nabla f_k\|\),即步长由半径决定的最速下降,两种框架此时本质相同。
  • 取 \(B_k=\nabla^2 f_k\):由于有信赖域约束,Hessian 非正定也无需特殊处理,子问题总有解(图 2.4)——信赖域牛顿法实践中非常有效(第 6 章)。
  • 取 \(B_k\) 为拟牛顿近似:信赖域拟牛顿法。

尺度(Scaling,PDF p.48–49)

  • 病态尺度(poorly scaled):\(x\) 在某方向的变化引起 \(f\) 的变化远大于另一方向。例:\(f(x)=10^9x_1^2+x_2^2\)。
  • 例:化学反应速率常数估计,粗略量级 \(x_1\approx10^{-10}\)、\(x_2\approx x_3\approx1\)、\(x_4\approx10^5\)。令 \(x=\mathrm{diag}(10^{-10},1,1,10^5)z\),以 \(z\) 为变量,最优 \(z\) 都在 1 的一个数量级内。这称为对角尺度化(diagonal scaling)。
  • 改变变量单位(米→毫米)就是在做(有时无意的)尺度变换。
  • 最速下降对尺度敏感,牛顿法不受影响(图 2.7:椭圆等高线拉长时最速下降方向几乎不降 \(f\);圆形时表现好;牛顿步两种情况都好,因为二次模型逼近好)。
  • 设计算法时应尽量在线搜索/信赖域策略和收敛检验中保持尺度不变性(scale invariance);一般线搜索比信赖域更容易保持。

收敛速度(Rates of Convergence,PDF p.49–50)

设 \(x_k\to x^*\):

  • Q-线性(Q-linear):存在 \(r\in(0,1)\),\(\dfrac{\|x_{k+1}-x^*\|}{\|x_k-x^*\|}\le r\) 对充分大 \(k\) 成立(2.21)。例:\(1+(0.5)^k\)。Q 代表 quotient(商)。
  • Q-超线性(Q-superlinear):\(\lim_{k\to\infty}\dfrac{\|x_{k+1}-x^*\|}{\|x_k-x^*\|}=0\)。例:\(1+k^{-k}\)。
  • Q-二次(Q-quadratic):\(\dfrac{\|x_{k+1}-x^*\|}{\|x_k-x^*\|^2}\le M\)(\(M>0\) 不必小于 1)。例:\(1+(0.5)^{2^k}\)。
  • 一般 Q-阶为 \(p>1\):\(\dfrac{\|x_{k+1}-x^*\|}{\|x_k-x^*\|^p}\le M\)。
  • 速度取决于 \(r\)(以及较弱地取决于 \(M\));无论常数如何,二次收敛序列最终总快于线性收敛序列。二次⇒超线性⇒线性。
  • 典型:拟牛顿法 Q-超线性;牛顿法 Q-二次;最速下降 Q-线性,病态时 \(r\) 接近 1。书中后文略去 “Q”。

R-收敛速度(R-Rates,PDF p.50–51)

  • R-线性(R 代表 root):存在非负序列 \(\{\nu_k\}\) 使 \(\|x_k-x^*\|\le\nu_k\) 且 \(\nu_k\) Q-线性收敛到 0(称误差被 \(\nu_k\) 控制,dominated)。例 (2.22):\(x_k=1+(0.5)^k\)(\(k\) 偶)、\(1\)(\(k\) 奇),前几项 \(2,1,1.25,1,1.03125,1,\dots\),R-线性收敛到 1。
  • 类似定义 R-超线性、R-二次。
  • 注意:R-收敛序列的误差可以在某些步增大(上例每隔一步误差增加),即使 R-阶任意高;Q-线性序列则要求充分大 \(k\) 后每步都减小。多数收敛分析关注 Q-收敛。

注释与参考(PDF p.51)

收敛速度的系统讨论见 Ortega & Rheinboldt。

习题概览(PDF p.51–53)

  • 2.1 Rosenbrock 函数 \(f=100(x_2-x_1^2)^2+(1-x_1)^2\)(2.23)的梯度、Hessian,证明 \(x^*=(1,1)\) 是唯一局部极小且 Hessian 正定。
  • 2.2 \(f=8x_1+12x_2+x_1^2-2x_2^2\) 唯一驻点是鞍点,画等高线。
  • 2.3 \(a^Tx\)、\(x^TAx\) 的梯度与 Hessian。
  • 2.4 \(\cos(1/x)\) 的二阶 Taylor 展开、\(\cos x\) 的三阶展开。
  • 2.5 \(f=\|x\|^2\),序列 \(x_k=(1+2^{-k})(\cos k,\sin k)^T\) 单调降但单位圆上每点都是极限点(说明“函数值下降≠收敛到一点”)。
  • 2.6 孤立局部极小必严格;2.7 凸函数全局极小点集为凸集。
  • 2.8 \(f=(x_1+x_2^2)^2\) 在 \((1,0)\) 沿 \((-1,1)\) 是下降方向,求 (2.9) 的极小。
  • 2.9 仿射变量变换 \(x=Sz+s\) 下 \(\nabla\tilde f=S^T\nabla f\),\(\nabla^2\tilde f=S^T\nabla^2 fS\);2.10 证明 SR1 与 BFGS 在适当 \(B_0\) 下尺度不变。
  • 2.11 病态尺度 ⇒ Hessian 病态。
  • 2.12–2.15 收敛速度判别:\(1/k\) 次线性;\(1+(0.5)^{2^k}\) Q-二次;\(1/k!\) 是否超线性/二次;交错序列的 Q/R 阶。

本章要点

  1. 局部/全局、严格、孤立极小的定义及相互关系(孤立⇒严格,反之不然)。
  2. Taylor 定理三种形式 (2.4)–(2.6) 是全书分析工具。
  3. 一阶必要(\(\nabla f=0\))、二阶必要(Hessian 半正定)、二阶充分(Hessian 正定⇒严格局部极小);凸函数驻点即全局极小。
  4. 两大框架:线搜索(先方向后步长)与信赖域(先半径后方向+步长)。
  5. 四类方向:最速下降、牛顿、拟牛顿(割线方程、SR1、BFGS 及其逆形式)、非线性 CG。
  6. 尺度问题:最速下降对尺度敏感,牛顿法尺度不变;对角尺度化。
  7. 收敛速度:Q-线性/超线性/二次及 R-收敛。

与量化交易的关联

  • 因子/模型参数估计:例 2.1 的非线性最小二乘与量化中拟合收益率曲线(Nelson–Siegel 模型 \(y(\tau)=\beta_0+\beta_1\frac{1-e^{-\tau/\lambda}}{\tau/\lambda}+\dots\))、波动率曲面(SVI 参数化)、GARCH 极大似然估计结构完全一致:参数少、数据多,每次目标求值要遍历全部样本。
  • 最优性条件用于检验求解结果:均值–方差无约束问题 \(\max_w\ \mu^Tw-\tfrac{\gamma}{2}w^T\Sigma w\) 的一阶条件 \(\mu=\gamma\Sigma w\) 给出 \(w^*=\gamma^{-1}\Sigma^{-1}\mu\);\(\Sigma\) 正定保证这是唯一全局解(定理 2.4/2.5)。若样本协方差奇异(股票数 > 样本期数),Hessian 只半正定,解不唯一,需收缩估计或正则化。
  • 非凸与局部极小:神经网络因子模型、多参数策略参数寻优的目标常多峰,要用多起点;本章说明局部算法无法保证全局解,回测中“参数调到最优”往往是在噪声上的局部极小。
  • 尺度问题:同时优化以“元”计的持仓和以“比例”计的权重、或因子暴露量级差异巨大时,梯度法收敛极慢;先标准化因子(z-score)或对变量做对角缩放,正是本章的对角尺度化。
  • 非光滑目标:\(\ell_1\) 交易成本 \(\sum|w_i-w_i^0|\)、最大回撤、CVaR 都是非光滑的,本章提示需要专门方法(通常转化为 LP 或用次梯度/近端方法)。
  • 收敛速度:实现中选择算法时,牛顿/拟牛顿在中小规模(几百个资产)上几十步收敛;最速下降在病态协方差矩阵(条件数常达 \(10^3\)–\(10^5\))下非常慢。

推荐习题

  • 2.1(Rosenbrock:熟练计算梯度/Hessian 并验证二阶充分条件,是后续数值实验的标准测试函数)。
  • 2.3(二次型梯度 \(\nabla(x^TAx)=2Ax\),组合优化推导必备)。
  • 2.5(理解函数值单调下降并不保证迭代点收敛)。
  • 2.9、2.10(仿射变换下的梯度/Hessian 变换与拟牛顿尺度不变性,对理解预条件和变量标准化很关键)。
  • 2.12–2.15(区分 Q-与 R-收敛阶)。

第 3 章 线搜索方法(Line Search Methods,PDF p.54–83)

3.0 引言(PDF p.55–56)

线搜索迭代:

\[x_{k+1}=x_k+\alpha_kp_k, \tag{3.1}\]
\(\alpha_k>0\) 为步长(step length)。成功取决于方向 \(p_k\) 和步长 \(\alpha_k\) 都选得好。多数方法要求 \(p_k\) 为下降方向(\(p_k^T\nabla f_k<0\)),且常取
\[p_k=-B_k^{-1}\nabla f_k, \tag{3.2}\]
\(B_k\) 对称非奇异:最速下降 \(B_k=I\);牛顿法 \(B_k=\nabla^2 f(x_k)\);拟牛顿法中 \(B_k\) 每步用低秩公式更新。若 \(B_k\) 正定,则 \(p_k^T\nabla f_k=-\nabla f_k^TB_k^{-1}\nabla f_k<0\),是下降方向。

3.1 步长(Step Length,PDF p.56–63)

权衡:希望 \(\alpha_k\) 使 \(f\) 大幅下降,又不愿花太多时间选择。理想选择是一元函数

\[\phi(\alpha)=f(x_k+\alpha p_k),\quad \alpha>0 \tag{3.3}\]
的全局极小点(图 3.1),但代价过高;即便求局部极小到中等精度也需要很多次函数/梯度求值。实用策略做非精确线搜索(inexact line search):以最小代价获得足够下降。典型过程分两阶段:包围阶段(bracketing)找到含合适步长的区间,二分或插值阶段(bisection/interpolation)在区间内算出好步长。

只要求 \(f\) 下降是不够的:图 3.2 中极小值 \(f^*=-1\),但若函数值序列为 \(\{5/k\}\),会收敛到 0 而非 \(-1\)。需要“充分下降”。

Wolfe 条件(PDF p.57–61)

充分下降条件(sufficient decrease,又称 Armijo 条件):

\[f(x_k+\alpha p_k)\le f(x_k)+c_1\alpha\nabla f_k^Tp_k,\quad c_1\in(0,1). \tag{3.4}\]
即下降量应与步长和方向导数成比例。右端记为线性函数 \(l(\alpha)\),斜率 \(c_1\nabla f_k^Tp_k<0\);由于 \(c_1<1\),对小正 \(\alpha\),\(l\) 位于 \(\phi\) 图像之上;条件为 \(\phi(\alpha)\le l(\alpha)\)(图 3.3)。实践中 \(c_1\) 很小,如 \(10^{-4}\)。

仅此不够:所有足够小的 \(\alpha\) 都满足它。为排除过短步长,引入曲率条件(curvature condition):

\[\nabla f(x_k+\alpha_kp_k)^Tp_k\ge c_2\nabla f_k^Tp_k,\quad c_2\in(c_1,1). \tag{3.5}\]
左端即 \(\phi'(\alpha_k)\),要求斜率大于 \(c_2\phi'(0)\)。直观:若斜率仍很负,说明继续前进还能显著下降;若斜率只略负或为正,说明此方向已难有更多下降,可终止(图 3.4)。典型值:牛顿/拟牛顿 \(c_2=0.9\);非线性共轭梯度 \(c_2=0.1\)。

Wolfe 条件(合称):

\[f(x_k+\alpha_kp_k)\le f(x_k)+c_1\alpha_k\nabla f_k^Tp_k, \tag{3.6a}\]
\[\nabla f(x_k+\alpha_kp_k)^Tp_k\ge c_2\nabla f_k^Tp_k, \tag{3.6b}\]
\(0<c_1<c_2<1\)。满足 Wolfe 条件的步长不一定接近 \(\phi\) 的极小点(图 3.5)。

强 Wolfe 条件(strong Wolfe):

\[f(x_k+\alpha_kp_k)\le f(x_k)+c_1\alpha_k\nabla f_k^Tp_k, \tag{3.7a}\]
\[|\nabla f(x_k+\alpha_kp_k)^Tp_k|\le c_2|\nabla f_k^Tp_k|. \tag{3.7b}\]
区别是不再允许 \(\phi'(\alpha_k)\) 太正,从而排除远离 \(\phi\) 驻点的点。

引理 3.1(Wolfe 步长存在性):设 \(f\) 连续可微,\(p_k\) 为 \(x_k\) 处下降方向,\(f\) 沿射线 \(\{x_k+\alpha p_k:\alpha>0\}\) 有下界。若 \(0<c_1<c_2<1\),则存在满足 Wolfe 条件 (3.6) 和强 Wolfe 条件 (3.7) 的步长区间。 证明:\(\phi\) 有下界而 \(l(\alpha)\) 斜率为负无下界,故 \(l\) 与 \(\phi\) 图像至少相交一次,设最小交点 \(\alpha'>0\),\(f(x_k+\alpha'p_k)=f(x_k)+\alpha'c_1\nabla f_k^Tp_k\)(3.8),小于 \(\alpha'\) 的步长都满足充分下降。由中值定理,存在 \(\alpha''\in(0,\alpha')\) 使 \(f(x_k+\alpha'p_k)-f(x_k)=\alpha'\nabla f(x_k+\alpha''p_k)^Tp_k\)(3.9)。合并得

\[\nabla f(x_k+\alpha''p_k)^Tp_k=c_1\nabla f_k^Tp_k>c_2\nabla f_k^Tp_k, \tag{3.10}\]
(因为 \(c_1<c_2\) 且 \(\nabla f_k^Tp_k<0\)),故 \(\alpha''\) 严格满足 Wolfe 两个条件,由光滑性其周围存在区间;又因 (3.10) 左端为负,强 Wolfe 也在该区间成立。

Wolfe 条件在广义上尺度不变:目标乘常数或做仿射变量变换不改变它们。可用于多数线搜索方法,对拟牛顿法实现尤其重要(第 8 章)。

Goldstein 条件(PDF p.61–62)

\[f(x_k)+(1-c)\alpha_k\nabla f_k^Tp_k\le f(x_k+\alpha_kp_k)\le f(x_k)+c\alpha_k\nabla f_k^Tp_k,\quad 0<c<\tfrac12. \tag{3.11}\]

右侧不等式为充分下降;左侧从下方控制步长(图 3.6)。缺点:左侧不等式可能排除 \(\phi\) 的所有极小点。与 Wolfe 收敛理论相似;常用于牛顿型方法,不适合需保持 Hessian 近似正定的拟牛顿法(因为 Goldstein 不保证 \(s_k^Ty_k>0\))。

充分下降与回溯(Sufficient Decrease and Backtracking,PDF p.61–63)

只用充分下降条件也可以,只要候选步长按回溯(backtracking)方式选取:

Procedure 3.1(回溯线搜索):

  1. 选 \(\bar\alpha>0\),\(\rho,c\in(0,1)\),置 \(\alpha\leftarrow\bar\alpha\);
  2. 当 \(f(x_k+\alpha p_k)>f(x_k)+c\alpha\nabla f_k^Tp_k\) 时,令 \(\alpha\leftarrow\rho\alpha\);
  3. 输出 \(\alpha_k=\alpha\)。
  • 牛顿与拟牛顿法中 \(\bar\alpha=1\);最速下降、CG 可取其他值。
  • 有限次试探必能终止(\(\alpha\) 足够小时充分下降成立)。
  • 实践中收缩因子 \(\rho\) 每次可变(如由带保护的插值确定),只需 \(\rho\in[\rho_{lo},\rho_{hi}]\),\(0<\rho_{lo}<\rho_{hi}<1\)。
  • 回溯保证:要么步长就是初始值 \(\bar\alpha\),要么足够短以满足充分下降但又不太短——因为被接受的 \(\alpha_k\) 与上一次被拒绝(太长)的 \(\alpha_k/\rho\) 只差一个因子。
  • 适合牛顿法(第 6 章),不太适合拟牛顿和 CG。

3.2 线搜索方法的收敛性(Convergence of Line Search Methods,PDF p.63–66)

全局收敛不仅需要好步长,也需要好方向。关键量是 \(p_k\) 与最速下降方向的夹角:

\[\cos\theta_k=\frac{-\nabla f_k^Tp_k}{\|\nabla f_k\|\|p_k\|}. \tag{3.12}\]

定理 3.2(Zoutendijk):考虑迭代 (3.1),\(p_k\) 为下降方向,\(\alpha_k\) 满足 Wolfe 条件 (3.6)。设 \(f\) 在 \(\mathbb{R}^n\) 有下界,在包含水平集 \(\mathcal{L}=\{x:f(x)\le f(x_0)\}\) 的开集 \(\mathcal{N}\) 上连续可微,且梯度 Lipschitz 连续:

\[\|\nabla f(x)-\nabla f(\tilde x)\|\le L\|x-\tilde x\|,\quad \forall x,\tilde x\in\mathcal{N}. \tag{3.13}\]
则
\[\sum_{k\ge0}\cos^2\theta_k\|\nabla f_k\|^2<\infty. \tag{3.14}\]
证明:由 (3.6b) 得 \((\nabla f_{k+1}-\nabla f_k)^Tp_k\ge(c_2-1)\nabla f_k^Tp_k\);由 Lipschitz 得 \((\nabla f_{k+1}-\nabla f_k)^Tp_k\le\alpha_kL\|p_k\|^2\)。合并得步长下界
\[\alpha_k\ge\frac{c_2-1}{L}\frac{\nabla f_k^Tp_k}{\|p_k\|^2}.\]
代入 (3.6a):\(f_{k+1}\le f_k-c_1\frac{1-c_2}{L}\frac{(\nabla f_k^Tp_k)^2}{\|p_k\|^2}=f_k-c\cos^2\theta_k\|\nabla f_k\|^2\),\(c=c_1(1-c_2)/L\)。累加得 \(f_{k+1}\le f_0-c\sum_{j=0}^k\cos^2\theta_j\|\nabla f_j\|^2\)(3.15),\(f\) 有下界故级数收敛。

  • 用 Goldstein 或强 Wolfe 条件也有类似结论。假设不苛刻:\(f\) 无下界则问题本身无定义;梯度 Lipschitz 常被局部收敛理论中的光滑性假设蕴含。
  • Zoutendijk 条件 (3.14) 蕴含 \(\cos^2\theta_k\|\nabla f_k\|^2\to0\)(3.16)。
  • 若方向与梯度夹角远离 \(90°\):\(\cos\theta_k\ge\delta>0\)(3.17),则
    \[\lim_{k\to\infty}\|\nabla f_k\|=0. \tag{3.18}\]
    特别地,用 Wolfe 或 Goldstein 线搜索的最速下降法,梯度收敛到零。
  • 全局收敛(globally convergent)在本书指满足 (3.18)。这是一般线搜索方法能得到的最强结果:只能保证被驻点吸引,不能保证收敛到极小点;需引入 Hessian 的负曲率信息才能加强为收敛到局部极小。
  • 牛顿型方法:若 \(B_k\) 正定且条件数一致有界 \(\|B_k\|\|B_k^{-1}\|\le M\),则 \(\cos\theta_k\ge1/M\)(3.19,习题 3.5),从而 \(\lim\|\nabla f_k\|=0\)(3.20)。即:\(B_k\) 正定、条件数有界、步长满足 Wolfe,则牛顿/拟牛顿法全局收敛。
  • 对 CG 等只能证明较弱的
    \[\liminf_{k\to\infty}\|\nabla f_k\|=0, \tag{3.21}\]
    即只有一个子列梯度趋于零。证明方法(反证):若 \(\|\nabla f_k\|\ge\gamma>0\)(3.22),则由 (3.16) 得 \(\cos\theta_k\to0\)(3.23);只需证明存在子列 \(\cos\theta_{k_j}\) 远离零即得矛盾。第 5 章用于非线性 CG。
  • 一般化:任何算法若 (i) 每步使 \(f\) 下降,(ii) 每 \(m\) 步做一次满足 Wolfe/Goldstein 的最速下降步,则 (3.20) 成立。其余 \(m-1\) 步可做“更好”的事;偶尔的最速下降步也许进展不大,但保证全局收敛。
  • 后续章节还会用到更强的“级数有界” (3.14),它迫使 \(\cos^2\theta_k\|\nabla f_k\|^2\) 足够快地趋于零。

3.3 收敛速度(Rate of Convergence,PDF p.66–75)

看似只要保证 \(p_k\) 不趋于与梯度正交(或定期走最速下降)即可。但若每步检查 \(\cos\theta_k\),小于 \(\delta\) 时把 \(p_k\) 拉向最速下降方向(角度检验),虽保证全局收敛却有两大弊端:(1) Hessian 病态时,好的方向本就可能几乎与梯度正交,\(\delta\) 选得不当会阻碍快速收敛;(2) 破坏拟牛顿法的不变性。快速收敛与全局收敛常常冲突:最速下降是全局收敛的典范但很慢;纯牛顿法在靠近解时极快,但远离解时步长可能不是下降方向。挑战在于兼得两者。

最速下降的收敛速度(PDF p.67–69)

理想情形:二次目标 + 精确线搜索。

\[f(x)=\tfrac12x^TQx-b^Tx,\quad Q\ \text{对称正定}. \tag{3.24}\]
\(\nabla f(x)=Qx-b\),\(x^*\) 为 \(Qx=b\) 的唯一解。沿 \(-\nabla f_k\) 的精确步长:
\[\alpha_k=\frac{\nabla f_k^T\nabla f_k}{\nabla f_k^TQ\nabla f_k}, \tag{3.25}\]
迭代 \(x_{k+1}=x_k-\dfrac{\nabla f_k^T\nabla f_k}{\nabla f_k^TQ\nabla f_k}\nabla f_k\)(3.26)。图 3.7:等高线为椭球,轴沿 \(Q\) 的特征向量,迭代呈之字形(zigzag)。

加权范数 \(\|x\|_Q^2=x^TQx\),利用 \(Qx^*=b\) 得

\[\tfrac12\|x-x^*\|_Q^2=f(x)-f(x^*), \tag{3.27}\]
并可推出(习题 3.7)
\[\|x_{k+1}-x^*\|_Q^2=\left\{1-\frac{(\nabla f_k^T\nabla f_k)^2}{(\nabla f_k^TQ\nabla f_k)(\nabla f_k^TQ^{-1}\nabla f_k)}\right\}\|x_k-x^*\|_Q^2. \tag{3.28}\]

定理 3.3:精确线搜索最速下降法用于强凸二次函数 (3.24) 时,

\[\|x_{k+1}-x^*\|_Q^2\le\left(\frac{\lambda_n-\lambda_1}{\lambda_n+\lambda_1}\right)^2\|x_k-x^*\|_Q^2, \tag{3.29}\]
\(0<\lambda_1\le\cdots\le\lambda_n\) 为 \(Q\) 的特征值(证明见 Luenberger;由 (3.28) 与 Kantorovich 不等式,习题 3.8)。

  • 函数值线性收敛。所有特征值相等(\(Q\) 为单位阵倍数,等高线为圆)时一步收敛。
  • 条件数 \(\kappa(Q)=\lambda_n/\lambda_1\) 越大,等高线越扁,之字形越严重,收敛越慢。(3.29) 虽为最坏界,但 \(n>2\) 时能准确反映实际行为。

定理 3.4(一般非线性):\(f\) 二阶连续可微,精确线搜索最速下降迭代收敛到 \(x^*\) 且 \(\nabla^2 f(x^*)\) 正定,则

\[f(x_{k+1})-f(x^*)\le\left(\frac{\lambda_n-\lambda_1}{\lambda_n+\lambda_1}\right)^2[f(x_k)-f(x^*)],\]
\(\lambda_i\) 为 \(\nabla^2f(x^*)\) 的特征值。非精确线搜索一般不会改善速度。 数值例:\(\kappa=800\),\(f(x_1)=1\),\(f(x^*)=0\),则 1000 次迭代后函数值仍约为 0.08。(验算:若按平方后的比率 \(r=(799/801)^2\approx0.99502\),\(r^{1000}\approx e^{-5}\approx0.0068\);原书的 0.08 对应未平方的 \((799/801)^{1000}\approx e^{-2.5}\approx0.082\)。编写教材时宜注明这一差异,结论“上千步仍未收敛到高精度”不变。)结论:即便 Hessian 条件尚可,最速下降也可能慢到不可接受。

拟牛顿法(Quasi-Newton Methods,PDF p.69–71)

方向 \(p_k=-B_k^{-1}\nabla f_k\)(3.30),\(B_k\) 对称正定,按拟牛顿公式更新;步长由满足 Wolfe/强 Wolfe 的非精确线搜索确定,且重要附加条件:线搜索总先试 \(\alpha=1\),满足 Wolfe 就接受(如 Procedure 3.1 中取 \(\bar\alpha=1\))。这一实现细节对快速收敛至关重要。

定理 3.5(Dennis–Moré):\(f\) 三阶连续可微,迭代 \(x_{k+1}=x_k+\alpha_kp_k\),\(p_k\) 为下降方向,\(\alpha_k\) 满足 Wolfe 条件且 \(c_1\le\frac12\)。若 \(x_k\to x^*\),\(\nabla f(x^*)=0\),\(\nabla^2 f(x^*)\) 正定,且

\[\lim_{k\to\infty}\frac{\|\nabla f_k+\nabla^2 f_kp_k\|}{\|p_k\|}=0, \tag{3.31}\]
则 (i) 存在 \(k_0\),\(k>k_0\) 时单位步长 \(\alpha_k=1\) 可被接受;(ii) 若 \(k>k_0\) 时都取 \(\alpha_k=1\),则 \(\{x_k\}\) 超线性收敛。

  • 若 \(c_1>\frac12\),线搜索会排除二次函数的极小点,单位步长可能不被接受。
  • 对拟牛顿方向,(3.31) 等价于
    \[\lim_{k\to\infty}\frac{\|(B_k-\nabla^2 f(x^*))p_k\|}{\|p_k\|}=0. \tag{3.32}\]
    惊喜结论:即使 \(B_k\) 不收敛到真 Hessian,只要沿搜索方向 \(p_k\) 越来越准,就能超线性收敛。(3.32) 是拟牛顿法超线性收敛的充要条件。

定理 3.6:\(f\) 三阶连续可微,迭代 \(x_{k+1}=x_k+p_k\)(步长恒为 1),\(p_k\) 由 (3.30) 给出,\(x_k\to x^*\),\(\nabla f(x^*)=0\),\(\nabla^2f(x^*)\) 正定。则 \(\{x_k\}\) 超线性收敛当且仅当 (3.32) 成立。 证明:先证 (3.32) ⇔ \(p_k-p_k^N=o(\|p_k\|)\)(3.33),其中 \(p_k^N=-\nabla^2f_k^{-1}\nabla f_k\): \(p_k-p_k^N=\nabla^2f_k^{-1}(\nabla^2f_kp_k+\nabla f_k)=\nabla^2f_k^{-1}(\nabla^2f_k-B_k)p_k=O(\|(\nabla^2f_k-B_k)p_k\|)=o(\|p_k\|)\)(用到 \(\|\nabla^2f_k^{-1}\|\) 在 \(x^*\) 附近有界);反向乘 \(\nabla^2f_k\) 即得。再用牛顿法二次收敛估计 (3.37): \(\|x_k+p_k-x^*\|\le\|x_k+p_k^N-x^*\|+\|p_k-p_k^N\|=O(\|x_k-x^*\|^2)+o(\|p_k\|)\),进而 \(\|p_k\|=O(\|x_k-x^*\|)\),得 \(\|x_k+p_k-x^*\|=o(\|x_k-x^*\|)\)。 第 8 章将说明拟牛顿法通常满足 (3.32)。

牛顿法(Newton's Method,PDF p.71–73)

\(p_k^N=-\nabla^2f_k^{-1}\nabla f_k\)(3.34)。Hessian 未必正定,故 \(p_k^N\) 未必下降,本章很多结论不适用。第 6 章给出两种全局化途径:线搜索(必要时修正 Hessian 使之正定)与信赖域。这里只讨论局部收敛:在 \(\nabla^2f(x^*)\) 正定的解附近,Hessian 也正定,若步长最终总为 1,则二次收敛。

定理 3.7:\(f\) 二阶可微,\(\nabla^2 f\) 在满足二阶充分条件的解 \(x^*\) 邻域内 Lipschitz 连续。迭代 \(x_{k+1}=x_k+p_k^N\)。则:(1) 初值足够接近 \(x^*\) 时迭代收敛到 \(x^*\);(2) 收敛速度二次;(3) \(\|\nabla f_k\|\) 二次收敛到零。 证明:

\[x_k+p_k^N-x^*=\nabla^2f_k^{-1}[\nabla^2f_k(x_k-x^*)-(\nabla f_k-\nabla f_*)]. \tag{3.35}\]
由 \(\nabla f_k-\nabla f_*=\int_0^1\nabla^2 f(x_k+t(x^*-x_k))(x_k-x^*)dt\),
\[\|\nabla^2f(x_k)(x_k-x^*)-(\nabla f_k-\nabla f(x^*))\|\le\int_0^1\|\nabla^2f(x_k)-\nabla^2f(x_k+t(x^*-x_k))\|\|x_k-x^*\|dt\le\tfrac12L\|x_k-x^*\|^2. \tag{3.36}\]
充分大 \(k\) 时 \(\|\nabla^2f_k^{-1}\|\le2\|\nabla^2f(x^*)^{-1}\|\),故
\[\|x_k+p_k^N-x^*\|\le L\|\nabla^2f(x^*)^{-1}\|\|x_k-x^*\|^2=\tilde L\|x_k-x^*\|^2. \tag{3.37}\]
归纳得收敛与二次速度。梯度:利用 \(\nabla f_k+\nabla^2f_kp_k^N=0\), \(\|\nabla f(x_{k+1})\|=\|\int_0^1\nabla^2f(x_k+tp_k^N)p_k^Ndt-\nabla^2f(x_k)p_k^N\|\le\tfrac12L\|p_k^N\|^2\le\tfrac12L\|\nabla^2f(x_k)^{-1}\|^2\|\nabla f_k\|^2\le2L\|\nabla^2f(x^*)^{-1}\|^2\|\nabla f_k\|^2\)。

牛顿方向使 (3.31) 中比值恒为零,故由定理 3.5,Wolfe(以及 Goldstein)条件对充分大的 \(k\) 接受单位步长。因此总先试单位步长的牛顿法实现最终 \(\alpha_k=1\),获得局部二次收敛。

坐标下降法(Coordinate Descent Methods,PDF p.73–75)

  • 依次以坐标方向 \(e_1,\dots,e_n\) 为搜索方向:固定其余变量,对一个变量极小化(或至少降低)\(f\),\(n\) 步后循环。又称交替变量法(method of alternating variables)。
  • 简单直观但实践中可能很低效(图 3.8:二维二次函数上,几步之后横竖移动都进展甚微)。
  • 带精确线搜索的坐标下降可无限迭代而永远不接近梯度为零的点(对比:最速下降保证 \(\|\nabla f_k\|\to0\))。推广(Powell):沿任意一组线性无关方向的循环搜索都不保证全局收敛。原因:梯度可能越来越垂直于坐标方向,\(\cos\theta_k\) 足够快趋零使 Zoutendijk 条件成立但 \(\nabla f_k\) 不趋零。
  • 即使收敛,通常比最速下降慢得多,变量越多差距越大。但仍有用:(1) 无需计算梯度;(2) 变量弱耦合(loosely coupled)时速度可接受。
  • 变体:“来回”顺序 \(e_1,\dots,e_{n-1},e_n,e_{n-1},\dots,e_1,e_2,\dots\);或一轮坐标步后沿首末点连线搜索(Hooke–Jeeves 等模式搜索方法基于此)。部分变体全局收敛。

3.4 步长选择算法(Step-Length Selection Algorithms,PDF p.75–81)

目标:极小化 \(\phi(\alpha)=f(x_k+\alpha p_k)\)(3.38)或只是找到满足 3.1 节终止条件的步长。假设 \(\phi'(0)<0\),只在 \(\alpha>0\) 上搜索。

  • 若 \(f\) 为凸二次 \(f=\tfrac12x^TQx+b^Tx+c\),沿射线的一维极小点有解析式:
    \[\alpha_k=-\frac{\nabla f_k^Tp_k}{p_k^TQp_k}. \tag{3.39}\]
  • 一般非线性函数需迭代。线搜索对所有非线性优化方法的稳健性和效率影响很大。
  • 只用函数值的算法可能低效(理论上要把极小点所在区间缩到很小);用梯度信息可判断是否已找到合适步长(Wolfe/Goldstein)。接近解时往往第一个试探步就满足条件。本节只讨论使用导数的算法。
  • 所有过程需初始估计 \(\alpha_0\),生成 \(\{\alpha_i\}\),终止于满足条件的步长或判定不存在。两阶段:包围阶段找区间 \([a,b]\);选择阶段(zoom)缩小区间并利用插值猜测极小点位置。
  • 记号:\(\alpha_k,\alpha_{k-1}\) 为外层迭代的步长;\(\alpha_i,\alpha_{i-1},\alpha_j\) 为线搜索内部试探步长;\(\alpha_0\) 为初始猜测。

插值(Interpolation,PDF p.76–78)

可视为 Procedure 3.1 的增强:生成递减序列 \(\alpha_i\),每个不比前一个小太多,直到满足充分下降

\[\phi(\alpha_k)\le\phi(0)+c_1\alpha_k\phi'(0). \tag{3.40}\]
\(c_1\) 很小(如 \(10^{-4}\)),此条件几乎只要求下降。设计目标是尽量少算导数。

  1. 若 \(\alpha_0\) 满足 (3.40),终止。否则 \([0,\alpha_0]\) 内含可接受步长。用 \(\phi(0)\)、\(\phi'(0)\)、\(\phi(\alpha_0)\) 构造二次插值
    \[\phi_q(\alpha)=\left(\frac{\phi(\alpha_0)-\phi(0)-\alpha_0\phi'(0)}{\alpha_0^2}\right)\alpha^2+\phi'(0)\alpha+\phi(0), \tag{3.41}\]
    其极小点
    \[\alpha_1=-\frac{\phi'(0)\alpha_0^2}{2[\phi(\alpha_0)-\phi(0)-\phi'(0)\alpha_0]}. \tag{3.42}\]
  2. 若 \(\alpha_1\) 仍不满足,用 \(\phi(0),\phi'(0),\phi(\alpha_0),\phi(\alpha_1)\) 构造三次插值 \(\phi_c(\alpha)=a\alpha^3+b\alpha^2+\alpha\phi'(0)+\phi(0)\),
    \[\begin{bmatrix}a\\b\end{bmatrix}=\frac{1}{\alpha_0^2\alpha_1^2(\alpha_1-\alpha_0)}\begin{bmatrix}\alpha_0^2&-\alpha_1^2\\-\alpha_0^3&\alpha_1^3\end{bmatrix}\begin{bmatrix}\phi(\alpha_1)-\phi(0)-\phi'(0)\alpha_1\\\phi(\alpha_0)-\phi(0)-\phi'(0)\alpha_0\end{bmatrix},\]
    极小点在 \([0,\alpha_1]\) 内:\(\alpha_2=\dfrac{-b+\sqrt{b^2-3a\phi'(0)}}{3a}\)。
  3. 必要时重复,用 \(\phi(0)\)、\(\phi'(0)\) 和最近两个 \(\phi\) 值做三次插值。
  4. 保护措施(safeguard):若 \(\alpha_i\) 离前一个太近或比前一个小太多,重置 \(\alpha_i=\alpha_{i-1}/2\),保证每步合理进展且最终 \(\alpha\) 不太小。

上述策略假设导数远比函数值贵。实际上常可几乎免费地与函数同时算方向导数(第 7 章),于是可用最近两个 \(\alpha\) 处的 \(\phi,\phi'\) 做三次插值(对曲率变化显著的函数是好模型)。设区间 \([a,b]\) 含合意步长,\(\alpha_{i-1},\alpha_i\) 在其中,插值 \(\phi(\alpha_{i-1}),\phi'(\alpha_{i-1}),\phi(\alpha_i),\phi'(\alpha_i)\) 的三次函数存在唯一;其在 \([a,b]\) 上的极小点要么在端点,要么在内部:

\[\alpha_{i+1}=\alpha_i-(\alpha_i-\alpha_{i-1})\left[\frac{\phi'(\alpha_i)+d_2-d_1}{\phi'(\alpha_i)-\phi'(\alpha_{i-1})+2d_2}\right], \tag{3.43}\]
\[d_1=\phi'(\alpha_{i-1})+\phi'(\alpha_i)-3\frac{\phi(\alpha_{i-1})-\phi(\alpha_i)}{\alpha_{i-1}-\alpha_i},\quad d_2=\left[d_1^2-\phi'(\alpha_{i-1})\phi'(\alpha_i)\right]^{1/2}.\]
迭代时丢弃 \(\alpha_{i-1}\) 或 \(\alpha_i\) 之一的数据,换入 \(\alpha_{i+1}\) 的数据;保留哪个取决于终止条件(见下文 Wolfe)。三次插值能使 (3.43) 二次收敛到极小步长。

初始步长(The Initial Step Length,PDF p.78)

  • 牛顿与拟牛顿:始终先试 \(\alpha_0=1\),以保证在满足条件时取单位步,发挥快速收敛性。
  • 最速下降、CG 等方向尺度不好的方法,需利用当前信息:
    • 假设本步一阶变化与上步相同:\(\alpha_0\nabla f_k^Tp_k=\alpha_{k-1}\nabla f_{k-1}^Tp_{k-1}\),即 \(\alpha_0=\alpha_{k-1}\dfrac{\nabla f_{k-1}^Tp_{k-1}}{\nabla f_k^Tp_k}\)。
    • 用 \(f(x_{k-1})\)、\(f(x_k)\)、\(\phi'(0)=\nabla f_k^Tp_k\) 做二次插值取极小:
      \[\alpha_0=\frac{2(f_k-f_{k-1})}{\phi'(0)}. \tag{3.44}\]
      若 \(x_k\to x^*\) 超线性,此比值趋于 1;取 \(\alpha_0\leftarrow\min(1,1.01\alpha_0)\),单位步最终总会被试并接受,从而在牛顿/拟牛顿中也能观察到超线性收敛。

满足 Wolfe 条件的线搜索算法(PDF p.78–81)

保证对任意 \(0<c_1<c_2<1\) 找到满足强 Wolfe (3.7) 的步长(假设 \(p\) 下降、\(f\) 沿 \(p\) 有下界)。两阶段:第一阶段从 \(\alpha_1\) 开始不断增大,直到找到可接受步长或包围区间;后者调用 zoom 不断缩小区间。\(\alpha_{\max}\) 为用户给定的最大步长。

Algorithm 3.2(线搜索算法):

α_0 ← 0;选 α_1 > 0 和 α_max;i ← 1
repeat
  计算 φ(α_i)
  if φ(α_i) > φ(0) + c1 α_i φ'(0)  或  [φ(α_i) ≥ φ(α_{i-1}) 且 i > 1]
      α* ← zoom(α_{i-1}, α_i);stop
  计算 φ'(α_i)
  if |φ'(α_i)| ≤ −c2 φ'(0)
      α* ← α_i;stop
  if φ'(α_i) ≥ 0
      α* ← zoom(α_i, α_{i-1});stop
  选 α_{i+1} ∈ (α_i, α_max);i ← i+1
end

试探步长单调增加,但传给 zoom 的参数顺序可能不同。依据:若 (i) \(\alpha_i\) 违反充分下降,或 (ii) \(\phi(\alpha_i)\ge\phi(\alpha_{i-1})\),或 (iii) \(\phi'(\alpha_i)\ge0\),则 \((\alpha_{i-1},\alpha_i)\) 内含满足强 Wolfe 的步长。外推选 \(\alpha_{i+1}\) 可用插值或取常数倍,关键是增长足够快以在有限步达到 \(\alpha_{\max}\)。

zoom(α_lo, α_hi) 的不变式:(a) \(\alpha_{lo}\) 与 \(\alpha_{hi}\) 之间含满足强 Wolfe 的步长;(b) \(\alpha_{lo}\) 是迄今满足充分下降的步长中函数值最小者;(c) \(\phi'(\alpha_{lo})(\alpha_{hi}-\alpha_{lo})<0\)。

Algorithm 3.3(zoom):

repeat
  用二次、三次插值或二分,在 α_lo 与 α_hi 之间求试探步长 α_j
  计算 φ(α_j)
  if φ(α_j) > φ(0) + c1 α_j φ'(0)  或  φ(α_j) ≥ φ(α_lo)
      α_hi ← α_j
  else
      计算 φ'(α_j)
      if |φ'(α_j)| ≤ −c2 φ'(0):  α* ← α_j;stop
      if φ'(α_j)(α_hi − α_lo) ≥ 0:  α_hi ← α_lo
      α_lo ← α_j
end

解释:\(\alpha_j\) 满足强 Wolfe 则终止;否则若它满足充分下降且函数值低于 \(\alpha_{lo}\),置 \(\alpha_{lo}\leftarrow\alpha_j\) 以维持 (b);若这破坏 (c),就把 \(\alpha_{hi}\) 设为旧 \(\alpha_{lo}\)。

实现要点:

  • 插值步要有保护,避免过于靠近区间端点;可利用插值多项式性质做有根据的猜测(Moré–Thuente 等)。
  • 接近解时 \(f(x_k)\) 与 \(f(x_{k-1})\) 可能在有限精度下无法区分,需加入停止测试:若若干次(典型 10 次)试探后无法降低函数值则停;或 \(x\) 的相对变化接近机器精度/用户阈值时停。
  • 完整实现很难写,建议用公开的高质量实现:Dennis–Schnabel、Lemaréchal、Fletcher,尤其 Moré–Thuente。
  • 强 Wolfe 与普通 Wolfe 代价:“松”线搜索(\(c_1=10^{-4}\),\(c_2=0.9\))两者工作量相近。强 Wolfe 的优势是减小 \(c_2\) 可直接控制搜索质量(迫使 \(\alpha\) 更接近局部极小),这对最速下降和非线性 CG 很重要,因此实现强 Wolfe 的例程适用面更广。

注释与参考(PDF p.81)

  • 终止条件详见 Ortega–Rheinboldt。Akaike 对二次函数上精确线搜索最速下降做了概率分析:\(n>2\) 时最坏界 (3.29) 对大多数初始点都成立;\(n=2\) 可闭式研究。
  • 一些线搜索方法(Goldfarb;Moré–Sorensen)在存在负曲率方向 \(p_-\)(\(p_-^T\nabla^2f(x_k)p_-<0\))时计算并组合它与最速下降方向,常做曲线回溯搜索,防止收敛到非极小驻点;由于难以确定两者相对权重,信赖域方法出现后此法不再流行。
  • 坐标下降收敛速度见 Luenberger。
  • 无导数线搜索:黄金分割(golden section)与 Fibonacci 搜索,存三个试探点确定含一维极小的区间,区别在于试探点生成方式。
  • 插值讨论遵循 Dennis–Schnabel;强 Wolfe 算法见 Fletcher。

习题概览(PDF p.82–83)

  • 3.1 编程:用回溯线搜索实现最速下降与牛顿法,极小化 Rosenbrock,\(\alpha_0=1\),打印每步步长;初值 \((1.2,1.2)\) 和较难的 \((-1.2,1)\)。
  • 3.2 若 \(0<c_2<c_1<1\),可能不存在 Wolfe 步长。
  • 3.3 证明 (3.39);3.4 \(c\le\frac12\) 时强凸二次的一维极小点满足 Goldstein。
  • 3.5 证明 \(\|Bx\|\ge\|x\|/\|B^{-1}\|\) 并推出 (3.19)。
  • 3.6 若 \(x_0-x^*\) 平行于 \(Q\) 的特征向量,最速下降一步到达解。
  • 3.7 推导 (3.28);3.8 Kantorovich 不等式 \(\dfrac{(x^Tx)^2}{(x^TQx)(x^TQ^{-1}x)}\ge\dfrac{4\lambda_n\lambda_1}{(\lambda_n+\lambda_1)^2}\),由此从 (3.28) 推出 (3.29)。
  • 3.9 编程:BFGS + 强 Wolfe 线搜索,验证 \(y_k^Ts_k>0\),解 Rosenbrock。
  • 3.10 证明 (3.41),并说明 \(\alpha_0\) 不满足充分下降时二次插值有正曲率且 \(\alpha_1<\frac{1}{2(1-c_1)}\)。
  • 3.11 构造使 \(\alpha_1\) 过小的例子;实践中取下界 \(\rho=0.1\):\(\alpha_i=\max(\rho\alpha_{i-1},\hat\alpha_i)\)。
  • 3.12 三次插值极小点在 \([0,\alpha_1]\);若 \(\phi(0)<\phi(\alpha_1)\) 则小于 \(\frac23\alpha_1\)。

本章要点

  1. 步长准则:Armijo 充分下降、曲率条件、(强) Wolfe、Goldstein;回溯法;Wolfe 步长存在性(引理 3.1)。
  2. Zoutendijk 定理:\(\sum\cos^2\theta_k\|\nabla f_k\|^2<\infty\);方向与梯度夹角远离 90° ⇒ 全局收敛(\(\|\nabla f_k\|\to0\));\(B_k\) 条件数有界 ⇒ 牛顿/拟牛顿全局收敛。
  3. 最速下降线性收敛,率 \(\left(\frac{\kappa-1}{\kappa+1}\right)^2\),病态时极慢、之字形。
  4. Dennis–Moré:\(\|(B_k-\nabla^2f^*)p_k\|/\|p_k\|\to0\) 是拟牛顿超线性收敛的充要条件,且此时单位步被 Wolfe 接受(需 \(c_1\le\frac12\)、总先试 \(\alpha=1\))。
  5. 牛顿法局部二次收敛(定理 3.7),梯度范数也二次收敛。
  6. 坐标下降可能不收敛、通常慢,但无需梯度、弱耦合时可用。
  7. 实用线搜索:二次/三次插值、保护、初始步长选取、Algorithm 3.2 + zoom(强 Wolfe),推荐 Moré–Thuente 实现。

与量化交易的关联

  • 系统实现:SciPy 的 minimize(method='BFGS'/'L-BFGS-B'/'CG') 内部就是 Wolfe 线搜索(scipy.optimize.line_search 即 Moré–Thuente 风格 + zoom);理解 \(c_1=10^{-4},c_2=0.9\)(拟牛顿)与 \(c_2=0.1\)(CG)的默认值,有助于调试“line search failed”报错(多因梯度实现错误或目标函数噪声导致无法下降)。
  • 风险模型病态与收敛:协方差矩阵条件数常很大(高相关资产、近似共线因子),用梯度下降拟合组合或风险平价时会呈之字形缓慢收敛,定理 3.3/3.4 给出定量解释;应改用牛顿/拟牛顿或预条件(第 5 章)。
  • 坐标下降:Lasso/Elastic Net 因子选择(glmnet)用的正是坐标下降,因为 \(\ell_1\) 项可分且坐标子问题有闭式(软阈值);本章指出普通光滑问题上坐标下降可能失败,但变量弱耦合时效果好——这解释了为什么在标准化、低相关的因子上 glmnet 很快,而高相关因子时变慢。风险平价的循环坐标算法(Griveau-Billion 等)也属此类。
  • 回测参数优化:回测收益关于参数常不光滑、带噪声,Wolfe 线搜索的梯度条件不可靠;本章说明无导数一维搜索(黄金分割)与 Hooke–Jeeves 模式搜索可作为替代。
  • 极大似然估计(GARCH、状态空间模型):用牛顿型方法时总先试单位步长,接近解后二次收敛——是计量软件快速收敛的原理。

推荐习题

  • 3.1、3.9(必做编程:回溯最速下降/牛顿与 BFGS+强 Wolfe 在 Rosenbrock 上的对比,直观感受收敛速度差异)。
  • 3.5(条件数有界 ⇒ 方向角有界)。
  • 3.7、3.8(Kantorovich 不等式与最速下降收敛率的推导)。
  • 3.2、3.4(Wolfe/Goldstein 参数约束的意义)。
  • 3.10、3.12(插值步长的性质,自己实现线搜索时有用)。

第 4 章 信赖域方法(Trust-Region Methods,PDF p.84–119)

4.0 引言与算法框架(PDF p.85–89)

  • 线搜索与信赖域都借助目标的二次模型,但用法不同:线搜索用模型产生方向再找步长;信赖域在当前点周围定义一个“相信模型能充分代表目标”的区域,在区域内取模型的近似极小作为步,同时选择方向和长度。步不可接受则缩小区域重求,方向一般随半径改变。
  • 区域大小至关重要:太小错失大步前进的机会;太大则模型极小点可能远离目标极小点。实践中依据之前迭代表现调整:模型可靠(步好、预测准)就逐步放大以允许更大胆的步;失败步说明模型在当前区域内不足,就缩小。
  • 图 4.1:当前点在弯曲山谷一端、极小点在另一端。基于同一模型的线搜索沿模型极小点方向搜索,即使最优步长也只有小幅下降;信赖域取虚线圆内的模型极小点,下降更显著。

模型(前两项与 Taylor 展开相同):

\[m_k(p)=f_k+\nabla f_k^Tp+\tfrac12p^TB_kp, \tag{4.1}\]
\(B_k\) 对称。由 (2.6),\(f(x_k+p)=f_k+\nabla f_k^Tp+\tfrac12p^T\nabla^2f(x_k+tp)p\)(4.2),模型误差 \(O(\|p\|^2)\);若 \(B_k=\nabla^2f(x_k)\),误差 \(O(\|p\|^3)\),称信赖域牛顿法(第 6 章)。本章只假设 \(B_k\) 对称且关于 \(k\) 一致有界。

信赖域子问题:

\[\min_{p\in\mathbb{R}^n}m_k(p)=f_k+\nabla f_k^Tp+\tfrac12p^TB_kp\quad\text{s.t. }\|p\|\le\Delta_k, \tag{4.3}\]
暂取欧氏范数,目标与约束(\(p^Tp\le\Delta_k^2\))都是二次的。若 \(B_k\) 正定且 \(\|B_k^{-1}\nabla f_k\|\le\Delta_k\),解就是无约束极小 \(p_k^B=-B_k^{-1}\nabla f_k\),称完全步(full step)。其他情况需计算,但只需近似解即可保证收敛与良好实践表现。

半径调整:定义比值

\[\rho_k=\frac{f(x_k)-f(x_k+p_k)}{m_k(0)-m_k(p_k)}, \tag{4.4}\]
分子为实际下降(actual reduction),分母为预测下降(predicted reduction,恒非负,因为 \(p=0\) 可行)。\(\rho_k<0\):目标上升,拒绝;\(\rho_k\approx1\):模型与函数吻合好,可扩大;\(\rho_k\) 为正但不接近 1:不变;接近 0 或为负:缩小。

Algorithm 4.1(信赖域):给定 \(\hat\Delta>0\)(步长总上界)、\(\Delta_0\in(0,\hat\Delta)\)、\(\eta\in[0,\frac14)\):

for k = 0,1,2,...
  (近似)求解 (4.3) 得 p_k;按 (4.4) 计算 ρ_k
  if ρ_k < 1/4:              Δ_{k+1} = (1/4)‖p_k‖
  else if ρ_k > 3/4 且 ‖p_k‖ = Δ_k:  Δ_{k+1} = min(2Δ_k, Δ̂)
  else:                       Δ_{k+1} = Δ_k
  if ρ_k > η:  x_{k+1} = x_k + p_k   else  x_{k+1} = x_k
end

只有当 \(\|p_k\|\) 真正到达边界时才放大半径;步严格在内部说明当前半径不妨碍进展,保持不变。

求解 (4.3) 的近似策略:先讲三种至少达到 Cauchy 点下降量的方法——dogleg(\(B_k\) 正定时适用)、二维子空间极小化(\(B_k\) 不定也可,需估计最负特征值)、Steihaug 方法(\(B_k\) 为大规模稀疏精确 Hessian 时最合适);再讲 Moré–Sorensen 的“近似精确”解法,基于解满足 \((B_k+\lambda I)p=-\nabla f_k\)(\(\lambda\ge0\)),寻找与半径相应的 \(\lambda\)。

4.1 Cauchy 点及相关算法(PDF p.89–97)

Cauchy 点(The Cauchy Point,PDF p.89–90)

与线搜索不需最优步长类似,信赖域全局收敛只需近似解位于信赖域内且使模型充分下降,后者用 Cauchy 点量化。

Algorithm 4.2(Cauchy 点计算):

  1. 求线性化子问题 \(p_k^S=\arg\min_p f_k+\nabla f_k^Tp\) s.t. \(\|p\|\le\Delta_k\)(4.5);
  2. 求 \(\tau_k=\arg\min_{\tau>0}m_k(\tau p_k^S)\) s.t. \(\|\tau p_k^S\|\le\Delta_k\)(4.6);
  3. \(p_k^C=\tau_kp_k^S\)。

闭式:\(p_k^S=-\dfrac{\Delta_k}{\|\nabla f_k\|}\nabla f_k\),

\[p_k^C=-\tau_k\frac{\Delta_k}{\|\nabla f_k\|}\nabla f_k, \tag{4.7}\]
\[\tau_k=\begin{cases}1,&\nabla f_k^TB_k\nabla f_k\le0;\\ \min\left(\dfrac{\|\nabla f_k\|^3}{\Delta_k\nabla f_k^TB_k\nabla f_k},1\right),&\text{否则}.\end{cases} \tag{4.8}\]
(曲率非正时沿负梯度模型单调降,走到边界;曲率为正时取一维凸二次的无约束极小或边界,先到者为准。)图 4.2 中 \(B_k\) 正定,Cauchy 点严格在区域内部。

Cauchy 步计算便宜(无需矩阵分解),是判断近似解是否可接受的关键:若每步使模型下降至少为 Cauchy 步下降的固定倍数,信赖域法即全局收敛。

改进 Cauchy 点(Improving on the Cauchy Point,PDF p.90–91)

总取 Cauchy 点等价于用特定步长的最速下降,表现差。Cauchy 点对 \(B_k\) 依赖弱(只用于确定步长);要想快速(如超线性)收敛,\(B_k\) 必须同时影响方向和长度。许多方法先算 Cauchy 点再改进,并设计成:当 \(B_k\) 正定且 \(\|p_k^B\|\le\Delta_k\) 时选完全步 \(p_k^B\)——当 \(B_k\) 为精确 Hessian 或拟牛顿近似时可期超线性收敛。

省略下标后的子问题:

\[\min_p m(p)\overset{\text{def}}{=}f+g^Tp+\tfrac12p^TBp\quad\text{s.t. }\|p\|\le\Delta, \tag{4.9}\]
解记为 \(p^*(\Delta)\)。

Dogleg 方法(PDF p.91–93)

  • \(B\) 正定时:\(\Delta\ge\|p^B\|\) 则 \(p^*(\Delta)=p^B\)(4.10);\(\Delta\) 很小时二次项影响小,\(p^*(\Delta)\approx-\Delta\dfrac{g}{\|g\|}\)(4.11);中间值时 \(p^*(\Delta)\) 沿一条曲线轨迹(图 4.3)。
  • Dogleg(狗腿):用两段折线代替曲线。第一段从原点到最速下降方向上的无约束极小点
    \[p^U=-\frac{g^Tg}{g^TBg}g, \tag{4.12}\]
    第二段从 \(p^U\) 到 \(p^B\):
    \[\tilde p(\tau)=\begin{cases}\tau p^U,&0\le\tau\le1,\\ p^U+(\tau-1)(p^B-p^U),&1\le\tau\le2.\end{cases} \tag{4.13}\]
    在信赖域约束下沿此路径极小化 \(m\)。

引理 4.1:\(B\) 正定时,(i) \(\|\tilde p(\tau)\|\) 关于 \(\tau\) 递增;(ii) \(m(\tilde p(\tau))\) 关于 \(\tau\) 递减。 证明(只需看 \(\tau\in[1,2]\)):(i) 令 \(h(\alpha)=\frac12\|p^U+\alpha(p^B-p^U)\|^2\), \(h'(\alpha)=-p^{U\,T}(p^U-p^B)+\alpha\|p^U-p^B\|^2\ge-p^{U\,T}(p^U-p^B)=\dfrac{g^Tg\,g^TB^{-1}g}{g^TBg}\left[1-\dfrac{(g^Tg)^2}{(g^TBg)(g^TB^{-1}g)}\right]\ge0\)(最后一步由 Cauchy–Schwarz,习题 4.6)。(ii) 令 \(\hat h(\alpha)=m(\tilde p(1+\alpha))\),\(\hat h'(\alpha)=(p^B-p^U)^T(g+Bp^U)+\alpha(p^B-p^U)^TB(p^B-p^U)\le(p^B-p^U)^T(g+Bp^B)=0\)。

推论:若 \(\|p^B\|\ge\Delta\),路径与边界恰交于一点,否则不相交。由于 \(m\) 沿路径递减:\(\|p^B\|\le\Delta\) 取 \(p^B\);否则取交点,解标量二次方程 \(\|p^U+(\tau-1)(p^B-p^U)\|^2=\Delta^2\) 得 \(\tau\)。无需搜索。 dogleg 可改造处理不定 \(B\),但意义不大(此时 \(p^B\) 不是 \(m\) 的无约束极小)。

二维子空间极小化(Two-Dimensional Subspace Minimization,PDF p.94)

  • \(B\) 正定时把搜索扩大到 \(p^U,p^B\) 张成(等价于 \(g,B^{-1}g\) 张成)的整个二维子空间:
    \[\min_p m(p)=f+g^Tp+\tfrac12p^TBp\quad\text{s.t. }\|p\|\le\Delta,\ p\in\mathrm{span}[g,B^{-1}g]. \tag{4.14}\]
    两个变量的问题,易解(习题 4.9)。Cauchy 点可行,故下降至少与 Cauchy 点一样多,全局收敛;整条 dogleg 路径也在此子空间内,故是 dogleg 的推广。
  • 优点:可直观、实用且理论上稳妥地处理不定 \(B\)(Byrd–Schnabel–Schultz):\(B\) 有负特征值时子空间改为
    \[\mathrm{span}[g,(B+\alpha I)^{-1}g],\quad \alpha\in(-\lambda_1,-2\lambda_1], \tag{4.15}\]
    \(\lambda_1\) 为 \(B\) 最负特征值(保证 \(B+\alpha I\) 正定;区间的灵活性允许用 Lanczos 等数值方法求 \(\alpha\))。若 \(\|(B+\alpha I)^{-1}g\|\le\Delta\),则放弃子空间搜索,取
    \[p=-(B+\alpha I)^{-1}g+v, \tag{4.16}\]
    \(v\) 满足 \(v^T(B+\alpha I)^{-1}g\le0\)(使 \(p\) 不往回缩,而继续大致沿 \(-(B+\alpha I)^{-1}g\) 方向走)。
  • \(B\) 有零特征值但无负特征值时,用 Cauchy 步。
  • 二维极小化得到的模型下降常接近精确解;主要计算量是一次 \(B\) 或 \(B+\alpha I\) 的分解,而近似精确方法通常需两三次分解。

Steihaug 方法(CG–Steihaug,PDF p.95–97)

前两种都需解一个含 \(B\) 的线性方程组,\(B\) 很大时昂贵。Steihaug 基于共轭梯度(CG),不需精确解线性方程组,又能改进 Cauchy 点。与标准 CG 的区别:在离开信赖域或遇到 \(B\) 的负曲率方向时终止。

Algorithm 4.3(CG–Steihaug):给定 \(\epsilon>0\),\(p_0=0\),\(r_0=g\),\(d_0=-r_0\);若 \(\|r_0\|<\epsilon\) 返回 \(p_0\)。

for j = 0,1,2,...
  if d_j^T B d_j ≤ 0:
     求 τ 使 p = p_j + τ d_j 在 (4.9) 中极小化 m 且 ‖p‖ = Δ;return p
  α_j = r_j^T r_j / d_j^T B d_j
  p_{j+1} = p_j + α_j d_j
  if ‖p_{j+1}‖ ≥ Δ:
     求 τ ≥ 0 使 p = p_j + τ d_j 满足 ‖p‖ = Δ;return p
  r_{j+1} = r_j + α_j B d_j
  if ‖r_{j+1}‖ < ε‖r_0‖: return p_{j+1}
  β_{j+1} = r_{j+1}^T r_{j+1} / r_j^T r_j
  d_{j+1} = r_{j+1} + β_{j+1} d_j      (PDF 文本此处为“+”,按标准 CG 及 d_0 = −r_0 的约定应为 d_{j+1} = −r_{j+1} + β_{j+1} d_j)
end
  • 对应关系:\(m\leftrightarrow\phi\),\(p\leftrightarrow x\),\(B\leftrightarrow A\),\(-g\leftrightarrow b\)。
  • 两个额外停止准则:遇到零/负曲率方向;\(p_{j+1}\) 越出信赖域。两种情况都把当前方向与边界的交点作为最终点。
  • 初始化 \(p_0=0\) 至关重要:第一步后 \(p_1=\alpha_0d_0=-\dfrac{g^Tg}{g^TBg}g\),恰为 Cauchy 点;CG 每步都降低 \(m\),故满足全局收敛的必要条件。
  • 另一关键性质:每个迭代点范数都比前一个大(也是 \(p_0=0\) 的结果),因此一到达边界就可以停止——之后不会再有信赖域内更低的点。

定理 4.2:Algorithm 4.3 生成的序列满足 \(0=\|p_0\|_2<\cdots<\|p_j\|_2<\|p_{j+1}\|_2<\cdots<\|p\|_2\le\Delta\)。 证明:先证 \(p_j^Tr_j=0\)(\(j\ge0\))和 \(p_j^Td_j>0\)(\(j\ge1\))。由 \(p_j=\sum_{i=0}^{j-1}\alpha_id_i\) 及 CG 的扩展子空间性质(\(r_j\perp d_i,i<j\),第 5 章)得 \(p_j^Tr_j=0\)。\(p_1^Td_1=(\alpha_0d_0)^T(r_1+\beta_1d_0)=\alpha_0\beta_1d_0^Td_0>0\)(4.17);归纳:\(p_{j+1}^Td_{j+1}=\beta_{j+1}p_{j+1}^Td_j=\beta_{j+1}(p_j^Td_j+\alpha_jd_j^Td_j)>0\)。于是 \(\|p_{j+1}\|^2=\|p_j\|^2+2\alpha_jp_j^Td_j+\alpha_j^2\|d_j\|^2>\|p_j\|^2\);若因负曲率或越界停止,则最终 \(\|p\|=\Delta\) 为最大可能长度。

解释:迭代点沿一条从 \(p_1\)(Cauchy 点)到 \(p\) 的插值路径,每步离起点越来越远;\(B\) 正定时可与 dogleg 路径类比(都从 \(p^C\) 走向 \(p^B\),直到碰到边界)。信赖域牛顿法取 \(B=\nabla^2f(x)\),迭代中可能不定,局部与全局收敛性都很好(第 6 章)。

4.2 子问题的近似精确解(Using Nearly Exact Solutions to the Subproblem,PDF p.97–107)

精确解的刻画(Characterizing Exact Solutions,PDF p.97–98)

上述方法不认真追求精确解,但利用了 \(B\) 的信息,成本低且全局收敛。\(n\) 不太大时,值得更充分利用模型:约用三次分解(dogleg/二维子空间只需一次)得到更好的近似。思路:解满足 \((B+\lambda I)p^*=-g\)(\(\lambda\ge0\)),对 \(\lambda\) 用一元牛顿法。

定理 4.3:\(p^*\) 是

\[\min_p m(p)=f+g^Tp+\tfrac12p^TBp\quad\text{s.t. }\|p\|\le\Delta \tag{4.18}\]
的全局解,当且仅当 \(p^*\) 可行且存在 \(\lambda\ge0\) 使
\[(B+\lambda I)p^*=-g, \tag{4.19a}\]
\[\lambda(\Delta-\|p^*\|)=0, \tag{4.19b}\]
\[B+\lambda I\ \text{半正定}. \tag{4.19c}\]

  • (4.19b) 为互补条件(complementarity):\(\lambda\) 与 \(\Delta-\|p^*\|\) 至少一个为零。解严格在内部(图 4.4 中 \(\Delta_1\))时 \(\lambda=0\),\(Bp^*=-g\) 且 \(B\) 半正定;在边界上(\(\Delta_2,\Delta_3\))\(\lambda\) 可为正。
  • 由 (4.19a),\(\lambda p^*=-Bp^*-g=-\nabla m(p^*)\):解与模型负梯度共线、垂直于模型等高线。
  • 注意:这是非凸二次在球上的全局最优性的充要条件——信赖域子问题是少数“非凸但可全局求解”的问题之一。

计算近似精确解(Calculating Nearly Exact Solutions,PDF p.98–102)

要么 \(\lambda=0\) 满足 (4.19a)(4.19c) 且 \(\|p\|\le\Delta\),要么定义 \(p(\lambda)=-(B+\lambda I)^{-1}g\)(\(\lambda\) 足够大使 \(B+\lambda I\) 正定,习题 4.10),求 \(\lambda>0\) 使

\[\|p(\lambda)\|=\Delta. \tag{4.20}\]
这是一维求根问题。设 \(B=Q\Lambda Q^T\),\(\Lambda=\mathrm{diag}(\lambda_1,\dots,\lambda_n)\),\(\lambda_1\le\cdots\le\lambda_n\),则对 \(\lambda\neq-\lambda_j\):
\[p(\lambda)=-Q(\Lambda+\lambda I)^{-1}Q^Tg=-\sum_{j=1}^n\frac{q_j^Tg}{\lambda_j+\lambda}q_j, \tag{4.21}\]
\[\|p(\lambda)\|^2=\sum_{j=1}^n\frac{(q_j^Tg)^2}{(\lambda_j+\lambda)^2}. \tag{4.22}\]

  • 在 \((-\lambda_1,\infty)\) 上 \(\|p(\lambda)\|\) 连续、非增,\(\lim_{\lambda\to\infty}\|p(\lambda)\|=0\)(4.23);当 \(q_j^Tg\ne0\) 时 \(\lim_{\lambda\to-\lambda_j}\|p(\lambda)\|=\infty\)(4.24)。因此(\(q_1^Tg\neq0\) 时)在 \((-\lambda_1,\infty)\) 上恰有一点 \(\lambda^*\) 使 \(\|p\|=\Delta\)(图 4.5)。
  • \(B\) 正定且 \(\|B^{-1}g\|\le\Delta\):\(\lambda=0\) 即解,无需搜索;\(B\) 正定但 \(\|B^{-1}g\|>\Delta\):在 \((0,\infty)\) 中搜索;\(B\) 不定且 \(q_1^Tg\neq0\):在 \((-\lambda_1,\infty)\) 中搜索。
  • 直接对 \(\phi_1(\lambda)=\|p(\lambda)\|-\Delta\)(4.25)用牛顿法不好:\(\lambda\) 略大于 \(-\lambda_1\) 时 \(\phi_1\approx\frac{C_1}{\lambda+\lambda_1}+C_2\),高度非线性,牛顿法不可靠或慢。改为
    \[\phi_2(\lambda)=\frac1\Delta-\frac{1}{\|p(\lambda)\|},\]
    在该区域 \(\phi_2\approx\frac1\Delta-\frac{\lambda+\lambda_1}{C_3}\),近似线性,牛顿法表现好(需保持 \(\lambda>-\lambda_1\),图 4.6)。牛顿迭代
    \[\lambda^{(\ell+1)}=\lambda^{(\ell)}-\frac{\phi_2(\lambda^{(\ell)})}{\phi_2'(\lambda^{(\ell)})}. \tag{4.26}\]

Algorithm 4.4(精确信赖域):给定 \(\lambda^{(0)}\)、\(\Delta>0\):

for ℓ = 0,1,2,...
  Cholesky 分解 B + λ^(ℓ) I = R^T R
  解 R^T R p_ℓ = −g,  R^T q_ℓ = p_ℓ
  λ^(ℓ+1) = λ^(ℓ) + (‖p_ℓ‖/‖q_ℓ‖)^2 · (‖p_ℓ‖ − Δ)/Δ        (4.27)
end

需加保护(如 \(\lambda^{(\ell)}<-\lambda_1\) 时 Cholesky 不存在)。每步主要工作是 Cholesky 分解;实用版不追求高精度,两三次迭代得到近似解即可。

困难情形(The Hard Case,PDF p.102–104)

  • 若最负特征值是重根(\(0>\lambda_1=\lambda_2=\cdots\)),只要对某个 \(\lambda_j=\lambda_1\) 的 \(j\) 有 \(q_j^Tg\neq0\),上述方法仍适用。
  • 若对所有 \(\lambda_j=\lambda_1\) 的 \(j\) 都有 \(q_j^Tg=0\),(4.24) 不成立,\((-\lambda_1,\infty)\) 内可能不存在 \(\|p(\lambda)\|=\Delta\) 的 \(\lambda\)(图 4.7),Moré–Sorensen 称之为困难情形(hard case)。由定理 4.3,\(\lambda\in[-\lambda_1,\infty)\),只能 \(\lambda=-\lambda_1\)。
  • 求 \(p\):仅删去 \(\lambda_j=\lambda_1\) 的项不够。\(B-\lambda_1I\) 奇异,取单位特征向量 \(z\)(\((B-\lambda_1I)z=0\),\(q_j^Tz=0\) 对 \(\lambda_j\neq\lambda_1\)),令
    \[p=\sum_{j:\lambda_j\ne\lambda_1}\frac{q_j^Tg}{\lambda_j+\lambda}q_j+\tau z, \tag{4.28}\]
    (PDF 文本此处无负号,与 (4.21) 不一致,应为 \(-\sum\frac{q_j^Tg}{\lambda_j+\lambda}q_j+\tau z\))则 \(\|p\|^2=\sum_{j:\lambda_j\neq\lambda_1}\frac{(q_j^Tg)^2}{(\lambda_j+\lambda)^2}+\tau^2\),总能选 \(\tau\) 使 \(\|p\|=\Delta\),且 (4.19) 在 \(\lambda=-\lambda_1\) 时成立。

定理 4.3 的证明(PDF p.104–107)

引理 4.4:\(m(p)=g^Tp+\frac12p^TBp\)(4.29),\(B\) 任意对称。(i) \(m\) 有极小当且仅当 \(B\) 半正定且 \(g\in\mathrm{range}(B)\);(ii) 极小唯一当且仅当 \(B\) 正定;(iii) \(B\) 半正定时,任何满足 \(Bp=-g\) 的 \(p\) 都是全局极小。 证明:(i) 充分性:存在 \(p\) 使 \(Bp=-g\),则 \(m(p+w)=m(p)+g^Tw+(Bp)^Tw+\frac12w^TBw=m(p)+\frac12w^TBw\ge m(p)\)。必要性:极小点处 \(\nabla m=Bp+g=0\) 故 \(g\in\mathrm{range}(B)\),\(\nabla^2m=B\) 半正定。(ii) 正定时 \(w\neq0\) 有 \(w^TBw>0\);反之若半正定非正定,存在 \(w\neq0\),\(Bw=0\),\(m(p+w)=m(p)\),不唯一。(iii) 由 (i) 的证明。 例:\(B=\mathrm{diag}(1,0,2)\)(特征值 0,1,2,奇异)。若 \(g_2=0\) 则 \(g\in\mathrm{range}(B)\),有极小;若 \(g_2\neq0\),沿 \(\alpha(0,-g_2,0)^T\),\(\alpha\uparrow\infty\) 可使 \(m\) 无限下降。

定理 4.3 证明:

  • 充分性:设存在 \(\lambda\ge0\) 满足 (4.19)。由引理 4.4(iii),\(p^*\) 是 \(\hat m(p)=g^Tp+\frac12p^T(B+\lambda I)p=m(p)+\frac\lambda2p^Tp\)(4.30)的全局极小,故 \(m(p)\ge m(p^*)+\frac\lambda2(p^{*T}p^*-p^Tp)\)(4.31)。由互补性 \(\lambda(\Delta^2-p^{*T}p^*)=0\),得 \(m(p)\ge m(p^*)+\frac\lambda2(\Delta^2-p^Tp)\ge m(p^*)\) 对一切 \(\|p\|\le\Delta\)。
  • 必要性:若 \(\|p^*\|<\Delta\),\(p^*\) 是 \(m\) 的无约束极小,\(\lambda=0\) 满足 (4.19)。若 \(\|p^*\|=\Delta\),\(p^*\) 也解 \(\min m(p)\) s.t. \(\|p\|=\Delta\);由约束优化一阶条件(见第 12 章 (12.30)),拉格朗日函数 \(\mathcal{L}(p,\lambda)=m(p)+\frac\lambda2(p^Tp-\Delta^2)\) 在 \(p^*\) 驻定,得 \((B+\lambda I)p^*=-g\)(4.32)。对所有 \(\|p\|=\Delta\),\(m(p)\ge m(p^*)+\frac\lambda2(p^{*T}p^*-p^Tp)\),代入 \(g\) 整理得 \(\frac12(p-p^*)^T(B+\lambda I)(p-p^*)\ge0\)(4.33);方向集 \(\{\pm(p-p^*)/\|p-p^*\|:\|p\|=\Delta\}\) 在单位球面上稠密,故 (4.19c) 成立。最后证 \(\lambda\ge0\):若只有负的 \(\lambda\) 满足 (4.19a)(4.19c),则由 (4.31),\(\|p\|\ge\|p^*\|=\Delta\) 时 \(m(p)\ge m(p^*)\),结合球内最优性,\(p^*\) 是 \(m\) 的无约束全局极小,由引理 4.4(i),\(Bp^*=-g\) 且 \(B\) 半正定,于是 \(\lambda=0\) 也满足条件,矛盾。

4.3 全局收敛(Global Convergence,PDF p.107–113)

Cauchy 点所得下降(PDF p.107–109)

全局收敛只需近似解达到 Cauchy 下降的固定比例。dogleg、二维子空间、Algorithm 4.3 的近似解满足

\[m_k(0)-m_k(p_k)\ge c_1\|\nabla f_k\|\min\left(\Delta_k,\frac{\|\nabla f_k\|}{\|B_k\|}\right),\quad c_1\in(0,1]. \tag{4.34}\]
“min”中二选一是信赖域方法的典型特征(源自半径约束);当 \(\Delta_k\) 为较小者时,类似 Wolfe 第一条件:模型下降与梯度和步长大小成比例。

引理 4.5:Cauchy 点满足 (4.34),\(c_1=\frac12\):

\[m_k(0)-m_k(p_k^C)\ge\tfrac12\|\nabla f_k\|\min\left(\Delta_k,\frac{\|\nabla f_k\|}{\|B_k\|}\right). \tag{4.35}\]
证明分三种情况:

  1. \(\nabla f_k^TB_k\nabla f_k\le0\):\(m_k(p_k^C)-m_k(0)=-\Delta_k\|\nabla f_k\|+\frac12\frac{\Delta_k^2}{\|\nabla f_k\|^2}\nabla f_k^TB_k\nabla f_k\le-\Delta_k\|\nabla f_k\|\)。
  2. 曲率为正且 \(\frac{\|\nabla f_k\|^3}{\Delta_k\nabla f_k^TB_k\nabla f_k}\le1\)(4.36):\(\tau\) 取内点,\(m_k(p_k^C)-m_k(0)=-\frac12\frac{\|\nabla f_k\|^4}{\nabla f_k^TB_k\nabla f_k}\le-\frac12\frac{\|\nabla f_k\|^4}{\|B_k\|\|\nabla f_k\|^2}=-\frac12\frac{\|\nabla f_k\|^2}{\|B_k\|}\)。
  3. 否则 \(\nabla f_k^TB_k\nabla f_k<\frac{\|\nabla f_k\|^3}{\Delta_k}\)(4.37),\(\tau=1\):\(m_k(p_k^C)-m_k(0)\le-\Delta_k\|\nabla f_k\|+\frac12\frac{\Delta_k^2}{\|\nabla f_k\|^2}\frac{\|\nabla f_k\|^3}{\Delta_k}=-\frac12\Delta_k\|\nabla f_k\|\)。

定理 4.6:若 \(\|p_k\|\le\Delta_k\) 且 \(m_k(0)-m_k(p_k)\ge c_2[m_k(0)-m_k(p_k^C)]\),则 \(p_k\) 满足 (4.34),\(c_1=c_2/2\)。特别地精确解满足 \(c_1=\frac12\)。dogleg、二维子空间、Algorithm 4.3 都有 \(m_k(p_k)\le m_k(p_k^C)\),故 \(c_1=\frac12\)。

收敛到驻点(Convergence to Stationary Points,PDF p.109–113)

两类结果:\(\eta=0\)(只要 \(f\) 下降就接受)时梯度有零极限点;\(\eta>0\)(实际下降至少为预测下降的小比例)时 \(\nabla f_k\to0\)。假设:\(B_k\) 一致有界;水平集 \(\{x:f(x)\le f(x_0)\}\)(4.38)有界;允许 \(\|p_k\|\le\gamma\Delta_k\),\(\gamma\ge1\)(4.39)。

定理 4.7(\(\eta=0\)):\(\|B_k\|\le\beta\),\(f\) 在水平集上连续可微且有下界,所有近似解满足 (4.34)(4.39),则

\[\liminf_{k\to\infty}\|\nabla f_k\|=0. \tag{4.40}\]
证明要点:

  • \(|\rho_k-1|=\left|\dfrac{m_k(p_k)-f(x_k+p_k)}{m_k(0)-m_k(p_k)}\right|\);由 Taylor,\(|m_k(p_k)-f(x_k+p_k)|\le(\beta/2)\|p_k\|^2+C_4(p_k)\|p_k\|\)(4.41),\(C_4\) 可通过限制 \(\|p_k\|\) 任意小。
  • 反设 \(\|\nabla f_k\|\ge\epsilon\)(\(k\ge K\))(4.42),则 \(m_k(0)-m_k(p_k)\ge c_1\epsilon\min(\Delta_k,\epsilon/\beta)\)(4.43),从而 \(|\rho_k-1|\le\dfrac{\gamma\Delta_k(\beta\gamma\Delta_k/2+C_4)}{c_1\epsilon\min(\Delta_k,\epsilon/\beta)}\)(4.44)。
  • 取 \(\bar\Delta\) 足够小使 \(\beta\gamma\Delta_k/2+C_4\le\frac{c_1\epsilon}{2\gamma}\)(4.45)且 \(\bar\Delta\le\epsilon/\beta\),则 \(\Delta_k\le\bar\Delta\) 时 \(|\rho_k-1|\le\frac14\),\(\rho_k>\frac34\),半径不减。故只有 \(\Delta_k\ge\bar\Delta\) 时才会缩小,\(\Delta_k\ge\min(\Delta_K,\bar\Delta/4)\)(4.46)。
  • 若有无穷子列 \(\rho_k\ge\frac14\),则 \(f(x_k)-f(x_{k+1})\ge\frac14c_1\epsilon\min(\Delta_k,\epsilon/\beta)\),\(f\) 有下界迫使该子列 \(\Delta_k\to0\),与 (4.46) 矛盾;否则最终 \(\rho_k<\frac14\),\(\Delta_k\) 每步缩为 \(\frac14\),也趋于 0,矛盾。

定理 4.8(\(\eta\in(0,\frac14)\)):\(\|B_k\|\le\beta\),\(f\) 在水平集上梯度 Lipschitz 连续且有下界,近似解满足 (4.34)(4.39),则

\[\lim_{k\to\infty}\nabla f_k=0. \tag{4.47}\]
证明要点(Schultz–Schnabel–Byrd):任取 \(\nabla f_m\neq0\),梯度 Lipschitz 常数 \(\beta_1\),令 \(\epsilon=\frac12\|\nabla f_m\|\),\(R=\epsilon/\beta_1\),球 \(\mathcal{B}(x_m,R)\) 内 \(\|\nabla f(x)\|\ge\epsilon\)。由定理 4.7 的推理,迭代点不能永远留在球内;设 \(x_{l+1}\) 是第一个离开者。对实际走步的迭代求和:若这些步 \(\Delta_k\le\epsilon/\beta\),则 \(f(x_m)-f(x_{l+1})\ge\eta c_1\epsilon\sum\Delta_k\ge\eta c_1\epsilon R=\eta c_1\epsilon^2/\beta_1\)(4.48);否则 \(\ge\eta c_1\epsilon^2/\beta\)(4.49)。\(f(x_k)\downarrow f^*>-\infty\)(4.50),于是
\[\|\nabla f_m\|^2\le\left[\tfrac14\eta c_1\min\left(\tfrac1\beta,\tfrac1{\beta_1}\right)\right]^{-1}(f(x_m)-f^*)\to0.\]

基于近似精确解的算法的收敛(PDF p.113–114)

Moré–Sorensen 的带保护求根牛顿法(含困难情形处理),终止准则保证

\[m(0)-m(p)\ge c_1(m(0)-m(p^*)),\quad \|p\|\le\gamma\Delta, \tag{4.51}\]
\(p^*\) 为精确解(实际不需知道 \(p^*\),由实用终止准则推出)。与 (4.34) 的主要区别:(4.51) 更好地利用二阶项 \(p^TBp\)。例:\(g=0\) 但 \(B\) 有负特征值(鞍点),(4.34) 右端为零(前述算法会停在鞍点),(4.51) 右端为正,迫使算法离开鞍点。只有当二阶项真实反映 \(f\) 时才值得如此关注——文献中只处理了 \(B=\nabla^2f(x)\) 的信赖域牛顿法。

定理 4.9:Algorithm 4.1 取 \(B_k=\nabla^2f(x_k)\),\(\eta\in(0,\frac14)\),近似解满足 (4.51),则 \(\lim\|\nabla f_k\|=0\)。若水平集紧,则要么在满足二阶必要条件的点终止,要么 \(\{x_k\}\) 在水平集中有满足二阶必要条件的极限点。(证明见 Moré–Sorensen §4。)——即信赖域牛顿法能避开鞍点,这是其相对线搜索方法的重要理论优势。

4.4 其他改进(Other Enhancements,PDF p.114–117)

尺度(Scaling,PDF p.114–116)

  • 病态尺度的拓扑症状:\(x^*\) 位于狭长山谷,附近等高线是高度偏心的椭圆。球形信赖域不合适:沿敏感方向模型只在短距离内可信,沿不敏感方向可信距离更长。信赖域形状应使边界上各点对模型的信心大致相同——用椭球信赖域:
    \[\|Dp\|\le\Delta, \tag{4.52}\]
    \(D\) 为正对角阵,尺度化子问题
    \[\min_p m_k(p)=f_k+\nabla f_k^Tp+\tfrac12p^TB_kp\quad\text{s.t. }\|Dp\|\le\Delta_k. \tag{4.53}\]
    \(f\) 对 \(x_i\) 敏感则 \(d_{ii}\) 大,反之接近零。\(D\) 可由二阶导 \(\partial^2f/\partial x_i^2\) 可靠地构造;可逐步变化,只要 \(d_{ii}\in[d_{lo},d_{hi}]\),\(0<d_{lo}\le d_{hi}<\infty\)。\(D\) 不必精确,无需复杂启发式。
  • 本章所有算法都可改为椭球信赖域,收敛理论只需表面修改。Algorithm 4.5(广义 Cauchy 点):\(p_k^S=\arg\min f_k+\nabla f_k^Tp\) s.t. \(\|Dp\|\le\Delta_k\)(4.54);\(\tau_k=\arg\min_{\tau>0}m_k(\tau p_k^S)\) s.t. \(\|\tau Dp_k^S\|\le\Delta_k\)(4.55);\(p_k^C=\tau_kp_k^S\)。闭式:
    \[p_k^S=-\frac{\Delta_k}{\|D^{-1}\nabla f_k\|}D^{-2}\nabla f_k, \tag{4.56}\]
    \[\tau_k=\begin{cases}1,&\nabla f_k^TD^{-2}B_kD^{-2}\nabla f_k\le0,\\ \min\left(\dfrac{\|D^{-1}\nabla f_k\|^3}{\Delta_k\nabla f_k^TD^{-2}B_kD^{-2}\nabla f_k},1\right),&\text{否则}.\end{cases} \tag{4.57}\]
  • 更简单:令 \(\tilde p=Dp\),子问题变为 \(\min_{\tilde p}f_k+\nabla f_k^TD^{-1}\tilde p+\frac12\tilde p^TD^{-1}B_kD^{-1}\tilde p\) s.t. \(\|\tilde p\|\le\Delta_k\),用 \(D^{-1}\nabla f_k\) 代替梯度、\(D^{-1}B_kD^{-1}\) 代替 \(B_k\) 即可套用球形理论与算法。

非欧氏信赖域(Non-Euclidean Trust Regions,PDF p.116–117)

可用 \(\|p\|_1\le\Delta_k\)、\(\|p\|_\infty\le\Delta_k\) 或其尺度版。无约束优化中无明显优势,但对约束问题有用。例:界约束问题 \(\min f(x)\) s.t. \(x\ge0\),子问题

\[\min_p m_k(p)\quad\text{s.t. }x_k+p\ge0,\ \|p\|\le\Delta_k. \tag{4.58}\]
欧氏范数下可行域是球与非负象限的交,几何上别扭;用 \(\infty\)-范数则是盒子 \(x_k+p\ge0\),\(-\Delta_ke\le p\le\Delta_ke\),可用标准二次规划技术求解。

注释与参考(PDF p.117)

Powell 首先证明类似定理 4.7 的结果(\(\eta=0\),假设更弱但分析更复杂);Moré 综述 1982 年前的发展,强调尺度化范数;近似精确解部分取自 Moré–Sorensen;Byrd–Schnabel–Schultz 给出非精确信赖域一般理论并提出二维子空间极小化,关注不定 \(B\) 的处理以获得比定理 4.7/4.8 更强的局部收敛;Dennis–Schnabel 综述。

习题概览(PDF p.117–119)

  • 4.1 \(f=10(x_2-x_1^2)^2+(1-x_1)^2\) 在 \((0,-1)\) 与 \((0,0.5)\) 处画二次模型等高线及 \(\Delta\) 从 0 到 2 的解族。
  • 4.2 编程:精确 Hessian 的 dogleg 法解 Rosenbrock,试验半径更新规则。
  • 4.3 编程:基于 CG–Steihaug 的信赖域牛顿法,解扩展 Rosenbrock \(\sum_{i=1}^n[(1-x_{2i-1})^2+10(x_{2i}-x_{2i-1}^2)^2]\),\(n=10,50\),报告每步是负曲率、触边界还是满足停止准则。
  • 4.4 迭代有界时存在梯度为零的极限点。
  • 4.5 验证 (4.8);4.6 Cauchy–Schwarz 推出 \(\gamma=\frac{\|g\|^4}{(g^TBg)(g^TB^{-1}g)}\le1\)。
  • 4.7 双狗腿(double-dogleg)路径:原点→\(p^C\)→\(\bar\gamma p^B\)(\(\bar\gamma\in(\gamma,1]\))→\(p^B\),证明范数单调增(后来测试显示与标准 dogleg 差别不大)。
  • 4.8 证明 (4.26) 与 (4.27) 等价;4.9 \(B\) 正定时二维子空间问题的解;4.10 任意对称 \(B\) 存在 \(\lambda\ge0\) 使 \(B+\lambda I\) 正定;4.11 验证 (4.56)(4.57)。
  • 4.12 反例:\(g=(-1/\epsilon,-1,-\epsilon^2)^T\),\(B=\mathrm{diag}(1/\epsilon^3,1,\epsilon^3)\),\(\Delta=0.5\),精确解模型下降 \(\frac38+O(\epsilon)\),二维子空间法仅 \(O(\epsilon)\)。

本章要点

  1. 信赖域子问题 (4.3);比值 \(\rho_k\) 驱动的半径更新(Algorithm 4.1)。
  2. Cauchy 点闭式 (4.7)(4.8) 及其下降估计 \(\frac12\|g\|\min(\Delta,\|g\|/\|B\|)\);达到 Cauchy 下降固定比例即全局收敛。
  3. 近似解法:dogleg(\(B\) 正定)、二维子空间(可处理不定)、CG–Steihaug(大规模,\(p_1\) 即 Cauchy 点,迭代范数单调增)。
  4. 精确解刻画(定理 4.3):\((B+\lambda I)p=-g\),互补,\(B+\lambda I\) 半正定;用 \(1/\|p(\lambda)\|\) 的牛顿法求 \(\lambda\)(Algorithm 4.4),困难情形加特征向量分量。
  5. 全局收敛:\(\eta=0\) 得 \(\liminf\|\nabla f_k\|=0\);\(\eta>0\) 得 \(\lim\nabla f_k=0\);精确 Hessian + 近似精确解能收敛到满足二阶必要条件的点(避开鞍点)。
  6. 椭球尺度化与 \(\infty\)-范数信赖域(界约束时得到盒约束 QP)。

与量化交易的关联

  • 非线性最小二乘与参数校准:Levenberg–Marquardt(原书第 10 章)本质是信赖域法,\((J^TJ+\lambda I)p=-J^Tr\) 正是定理 4.3 的形式;期权定价模型校准(Heston、SABR 拟合隐含波动率曲面)普遍用 LM 或 scipy.optimize.least_squares(method='trf')(信赖域反射法),理解 \(\lambda\) 与半径的对应有助于调参。
  • 非凸目标的稳健性:风险平价、最大分散化、含非线性冲击成本的执行优化等可能 Hessian 不定,信赖域法无需修正 Hessian 即可处理负曲率,且能避开鞍点(定理 4.9)。SciPy 中 trust-ncg、trust-exact、trust-krylov 分别对应 Steihaug、Moré–Sorensen 精确法与 Lanczos 版本。
  • 岭回归的联系:\((B+\lambda I)p=-g\) 与岭回归 \((X^TX+\lambda I)\beta=X^Ty\) 同构——“\(\ell_2\) 范数约束 \(\|\beta\|\le\Delta\)”与“岭惩罚”通过拉格朗日乘子 \(\lambda\) 一一对应(图 4.5 的单调关系),这给因子收益回归中的收缩估计提供了几何解释。
  • 界约束与盒形信赖域:组合权重上下界 \(l\le w\le u\) 用 \(\infty\)-范数信赖域时子问题变成盒约束 QP,这就是 L-BFGS-B/TRF 类算法处理权重界的基本思路。
  • 尺度:不同资产/因子量级差异大时用椭球信赖域(对角缩放),类似先对变量标准化。

推荐习题

  • 4.2、4.3(实现 dogleg 与 CG–Steihaug,体会信赖域在 Rosenbrock 与大规模问题上的行为,是实现自研优化器的基础)。
  • 4.6、4.7(Cauchy–Schwarz 与 dogleg 路径单调性)。
  • 4.8(推导精确信赖域的 \(\lambda\) 牛顿迭代)。
  • 4.10、4.12(理解 \(B+\lambda I\) 正定化与二维子空间法的局限)。

第 5 章 共轭梯度法(Conjugate Gradient Methods,PDF p.120–153)

5.0 引言(PDF p.121–122)

CG 有双重用途:求解大型线性方程组的最有用技术之一;可推广到非线性优化。分别称线性 CG 和非线性 CG。

  • 线性 CG 由 Hestenes 与 Stiefel 于 1950 年代提出,作为求解正定系数矩阵线性方程组的迭代法,是高斯消元的替代,非常适合大规模问题。其性能与系数矩阵特征值分布紧密相关;通过预条件(preconditioning)改善分布可显著加速,预条件是实用 CG 的关键。
  • 首个非线性 CG 由 Fletcher 与 Reeves 于 1960 年代提出,是最早的大规模非线性优化技术之一;关键特征:不需存储矩阵,比最速下降快。

5.1 线性共轭梯度法(The Linear Conjugate Gradient Method,PDF p.122–140)

求解

\[Ax=b,\quad A\ \text{对称正定}, \tag{5.1}\]
等价于极小化
\[\phi(x)=\tfrac12x^TAx-b^Tx, \tag{5.2}\]
二者同解且唯一。\(\phi\) 的梯度等于方程组的残差:
\[\nabla\phi(x)=Ax-b\overset{\text{def}}{=}r(x). \tag{5.3}\]

共轭方向法(Conjugate Direction Methods,PDF p.122–127)

共轭(conjugacy):非零向量组 \(\{p_0,\dots,p_l\}\) 关于对称正定 \(A\) 共轭,若

\[p_i^TAp_j=0,\quad\forall i\neq j. \tag{5.4}\]
共轭向量组必线性无关(习题 5.2),故 \(A\) 至多有 \(n\) 个共轭方向。

共轭方向法:给定 \(x_0\) 和共轭方向组 \(\{p_0,\dots,p_{n-1}\}\),

\[x_{k+1}=x_k+\alpha_kp_k, \tag{5.5}\]
\[\alpha_k=-\frac{r_k^Tp_k}{p_k^TAp_k}\quad(\text{沿 }p_k\text{ 的精确一维极小,见 (3.39)}). \tag{5.6}\]

定理 5.1:对任意 \(x_0\),共轭方向法至多 \(n\) 步收敛到 (5.1) 的解 \(x^*\)。 证明:方向线性无关张成 \(\mathbb{R}^n\),\(x^*-x_0=\sum\sigma_kp_k\),左乘 \(p_k^TA\) 由共轭性得 \(\sigma_k=\dfrac{p_k^TA(x^*-x_0)}{p_k^TAp_k}\)(5.7)。又 \(x_k=x_0+\sum_{i<k}\alpha_ip_i\),故 \(p_k^TA(x_k-x_0)=0\),从而 \(p_k^TA(x^*-x_0)=p_k^TA(x^*-x_k)=p_k^T(b-Ax_k)=-p_k^Tr_k\),即 \(\sigma_k=\alpha_k\)。

几何解释:若 \(A\) 对角,等高线椭圆轴与坐标轴对齐,依次沿 \(e_1,\dots,e_n\) 一维极小化 \(n\) 步即得解(图 5.1);\(A\) 非对角时坐标轮换不再有限步收敛(图 5.2)。变换 \(\hat x=S^{-1}x\)(5.8),\(S=[p_0\ p_1\cdots p_{n-1}]\),则 \(\hat\phi(\hat x)=\phi(S\hat x)=\frac12\hat x^T(S^TAS)\hat x-(S^Tb)^T\hat x\),由共轭性 \(S^TAS\) 对角,在 \(\hat x\) 空间做坐标搜索等价于在 \(x\) 空间沿 \(p_i\) 搜索——共轭方向法就是“在使 Hessian 对角化的坐标系中做坐标下降”。

另一性质:Hessian 对角时每次坐标极小化正确确定解的一个分量,即 \(k\) 步后已在 \(e_1,\dots,e_k\) 张成的子空间上极小化。一般情形由下面定理给出。利用

\[r_{k+1}=r_k+\alpha_kAp_k. \tag{5.9}\]

定理 5.2(扩展子空间极小化,Expanding Subspace Minimization):对任意 \(x_0\),共轭方向法生成的 \(x_k\) 满足

\[r_k^Tp_i=0,\quad i=0,\dots,k-1, \tag{5.10}\]
且 \(x_k\) 是 \(\phi\) 在仿射集 \(\{x:x=x_0+\mathrm{span}\{p_0,\dots,p_{k-1}\}\}\)(5.11)上的极小点。 证明:\(\tilde x\) 在 (5.11) 上极小 ⇔ \(r(\tilde x)^Tp_i=0\),\(i<k\)(\(h(\sigma)=\phi(x_0+\sum\sigma_ip_i)\) 严格凸二次,令偏导为零并用链式法则)。归纳:\(r_1^Tp_0=0\) 由 \(\alpha_0\) 是一维极小;设 \(r_{k-1}^Tp_i=0\)(\(i\le k-2\)),由 (5.9),\(p_{k-1}^Tr_k=p_{k-1}^Tr_{k-1}+\alpha_{k-1}p_{k-1}^TAp_{k-1}=0\)(由 \(\alpha_{k-1}\) 定义),\(p_i^Tr_k=p_i^Tr_{k-1}+\alpha_{k-1}p_i^TAp_{k-1}=0\)(归纳假设 + 共轭)。

“当前残差与所有之前的搜索方向正交” (5.10) 在本章被大量使用。

共轭方向组的选取:\(A\) 的特征向量既正交又共轭,但大规模问题求全部特征向量太贵;修改 Gram–Schmidt 正交化可生成共轭方向,但要存储整个方向组,也贵。

CG 的基本性质(Basic Properties,PDF p.127–130)

CG 是一种特殊的共轭方向法:生成新方向 \(p_k\) 只需前一个方向 \(p_{k-1}\),就自动与所有更早方向共轭——存储和计算都很少。

\[p_k=-r_k+\beta_kp_{k-1}, \tag{5.12}\]
由 \(p_{k-1}^TAp_k=0\) 得 \(\beta_k=\dfrac{r_k^TAp_{k-1}}{p_{k-1}^TAp_{k-1}}\);\(p_0\) 取最速下降方向。

Algorithm 5.1(CG 初步版):\(r_0=Ax_0-b\),\(p_0=-r_0\);当 \(r_k\neq0\): \(\alpha_k=-\frac{r_k^Tp_k}{p_k^TAp_k}\)(5.13a);\(x_{k+1}=x_k+\alpha_kp_k\)(5.13b);\(r_{k+1}=Ax_{k+1}-b\)(5.13c);\(\beta_{k+1}=\frac{r_{k+1}^TAp_k}{p_k^TAp_k}\)(5.13d);\(p_{k+1}=-r_{k+1}+\beta_{k+1}p_k\)(5.13e)。

Krylov 子空间(Krylov subspace):

\[\mathcal{K}(r_0;k)\overset{\text{def}}{=}\mathrm{span}\{r_0,Ar_0,\dots,A^kr_0\}. \tag{5.14}\]

定理 5.3:若第 \(k\) 个迭代点不是解,则

\[r_k^Tr_i=0,\ i=0,\dots,k-1; \tag{5.15}\]
\[\mathrm{span}\{r_0,\dots,r_k\}=\mathrm{span}\{r_0,Ar_0,\dots,A^kr_0\}; \tag{5.16}\]
\[\mathrm{span}\{p_0,\dots,p_k\}=\mathrm{span}\{r_0,Ar_0,\dots,A^kr_0\}; \tag{5.17}\]
\[p_k^TAp_i=0,\ i=0,\dots,k-1. \tag{5.18}\]
因此 \(\{x_k\}\) 至多 \(n\) 步收敛到 \(x^*\)。 证明(归纳):

  • (5.16):由归纳假设 \(r_k,p_k\in\mathcal{K}(r_0;k)\),故 \(Ap_k\in\mathrm{span}\{Ar_0,\dots,A^{k+1}r_0\}\)(5.19),由 (5.9) 得 \(r_{k+1}\in\mathcal{K}(r_0;k+1)\),“⊂”成立;反向:\(A^{k+1}r_0=A(A^kr_0)\in\mathrm{span}\{Ap_0,\dots,Ap_k\}\),而 \(Ap_i=(r_{i+1}-r_i)/\alpha_i\),故 \(A^{k+1}r_0\in\mathrm{span}\{r_0,\dots,r_{k+1}\}\)。
  • (5.17):\(\mathrm{span}\{p_0,\dots,p_{k+1}\}=\mathrm{span}\{p_0,\dots,p_k,r_{k+1}\}\)(由 5.13e)\(=\mathrm{span}\{r_0,\dots,A^kr_0,r_{k+1}\}=\mathrm{span}\{r_0,\dots,r_{k+1}\}=\mathcal{K}(r_0;k+1)\)。
  • (5.18):\(p_{k+1}^TAp_i=-r_{k+1}^TAp_i+\beta_{k+1}p_k^TAp_i\)(5.20)。\(i=k\) 时由 \(\beta\) 的定义为零。\(i\le k-1\):由定理 5.2,\(r_{k+1}^Tp_i=0\)(\(i\le k\))(5.21);又 \(Ap_i\in\mathrm{span}\{Ar_0,\dots,A^{i+1}r_0\}\subset\mathrm{span}\{p_0,\dots,p_{i+1}\}\)(5.22),故 \(r_{k+1}^TAp_i=0\);第二项由归纳假设为零。
  • (5.15)(非归纳):由 (5.10) \(r_k^Tp_i=0\),而 \(p_i=-r_i+\beta_ip_{i-1}\) 推出 \(r_i\in\mathrm{span}\{p_i,p_{i-1}\}\),故 \(r_k^Tr_i=0\)。

注意:证明依赖 \(p_0=-r_0\),换别的 \(p_0\) 结论不成立。由于是梯度(残差)相互正交、搜索方向关于 \(A\) 共轭,“共轭梯度”其实是个误称。

实用形式(A Practical Form,PDF p.131–132)

利用 (5.13e)(5.10):\(\alpha_k=\dfrac{r_k^Tr_k}{p_k^TAp_k}\);利用 \(\alpha_kAp_k=r_{k+1}-r_k\) 及正交性:\(\beta_{k+1}=\dfrac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}\)。

Algorithm 5.2(CG 标准形式):

给定 x_0;r_0 ← A x_0 − b;p_0 ← −r_0;k ← 0
while r_k ≠ 0
  α_k ← r_kᵀr_k / p_kᵀA p_k              (5.23a)
  x_{k+1} ← x_k + α_k p_k                  (5.23b)
  r_{k+1} ← r_k + α_k A p_k                (5.23c)
  β_{k+1} ← r_{k+1}ᵀr_{k+1} / r_kᵀr_k      (5.23d)
  p_{k+1} ← −r_{k+1} + β_{k+1} p_k         (5.23e)
  k ← k+1
end
  • 只需保存最近两步的 \(x,r,p\),可覆盖旧值。每步主要计算:一次矩阵–向量乘 \(Ap_k\)、两个内积、三个向量加(\(O(n)\) 浮点运算);矩阵–向量乘的代价依问题而定。
  • 只推荐用于大规模问题;小问题宜用高斯消元或 SVD 等分解法(对舍入误差更不敏感)。大规模时 CG 的优点:不改变系数矩阵,不产生填充(fill-in);且有时收敛非常快。

收敛速度(Rate of Convergence,PDF p.132–137)

精确算术下至多 \(n\) 步终止;更重要的是特征值分布有利时远少于 \(n\) 步。

由 (5.23b)(5.17):

\[x_{k+1}=x_0+\gamma_0r_0+\gamma_1Ar_0+\cdots+\gamma_kA^kr_0=x_0+P_k^*(A)r_0, \tag{5.24–5.25}\]
\(P_k^*\) 为 \(k\) 次多项式。用 \(A\)-范数 \(\|z\|_A^2=z^TAz\)(5.26),有 \(\frac12\|x-x^*\|_A^2=\phi(x)-\phi(x^*)\)(5.27)。由定理 5.2,\(P_k^*\) 解
\[\min_{P_k}\|x_0+P_k(A)r_0-x^*\|_A, \tag{5.28}\]
即在所有前 \(k\) 步限于 Krylov 子空间的方法中,CG 使 \(A\)-范数误差最小(最优性)。

由 \(r_0=A(x_0-x^*)\),\(x_{k+1}-x^*=[I+P_k^*(A)A](x_0-x^*)\)(5.29)。设 \(A\) 特征值 \(0<\lambda_1\le\cdots\le\lambda_n\)、正交特征向量 \(v_i\),\(x_0-x^*=\sum\xi_iv_i\)(5.30),\(P_k(A)v_i=P_k(\lambda_i)v_i\),故

\[\|x_{k+1}-x^*\|_A^2=\sum_{i=1}^n\lambda_i[1+\lambda_iP_k^*(\lambda_i)]^2\xi_i^2, \tag{5.31}\]
\[\|x_{k+1}-x^*\|_A^2\le\min_{P_k}\max_{1\le i\le n}[1+\lambda_iP_k(\lambda_i)]^2\,\|x_0-x^*\|_A^2. \tag{5.32}\]
收敛速度由 \(\min_{P_k}\max_i[1+\lambda_iP_k(\lambda_i)]^2\)(5.33)刻画:找一个在 0 处取 1、在所有特征值处尽量小的 \(k+1\) 次多项式。

定理 5.4:若 \(A\) 只有 \(r\) 个不同特征值,则 CG 至多 \(r\) 步终止于解。 证明:设不同特征值 \(\tau_1<\cdots<\tau_r\),令 \(Q_r(\lambda)=\dfrac{(-1)^r}{\tau_1\cdots\tau_r}(\lambda-\tau_1)\cdots(\lambda-\tau_r)\),则 \(Q_r(\lambda_i)=0\),\(Q_r(0)=1\);\(\bar P_{r-1}(\lambda)=(Q_r(\lambda)-1)/\lambda\) 为 \(r-1\) 次多项式,代入 (5.33)(\(k=r-1\))得常数为 0,故 \(x_r=x^*\)。

定理 5.5(Luenberger):

\[\|x_{k+1}-x^*\|_A^2\le\left(\frac{\lambda_{n-k}-\lambda_1}{\lambda_{n-k}+\lambda_1}\right)^2\|x_0-x^*\|_A^2. \tag{5.34}\]
思路:取 \(Q_{k+1}(\lambda)=1+\lambda\bar P_k(\lambda)\),使其根在最大的 \(k\) 个特征值 \(\lambda_n,\dots,\lambda_{n-k+1}\) 以及 \(\lambda_1\) 与 \(\lambda_{n-k}\) 的中点;其在其余特征值上的最大值恰为 \((\lambda_{n-k}-\lambda_1)/(\lambda_{n-k}+\lambda_1)\)。

应用(特征值聚类):若 \(A\) 有 \(m\) 个大特征值,其余 \(n-m\) 个聚集在 1 附近(图 5.3),令 \(\epsilon=\lambda_{n-m}-\lambda_1\),则 \(m+1\) 步后 \(\|x_{m+1}-x^*\|\approx\epsilon\|x_0-x^*\|_A\),\(\epsilon\) 小时 \(m+1\) 步就得到好的近似。

  • 图 5.4 数值例:5 个大特征值,其余在 \([0.95,1.05]\);定理 5.5 预测第 6 步误差陡降,实际第 5 步就降了(定理只是上界);第 7 步再次陡降(几乎只有 6 个不同特征值,因 1 附近略分散而晚一步)。对比特征值随机均匀分布的问题,收敛更慢且更均匀。
  • 一般结论:特征值成 \(r\) 个簇时,CG 约 \(r\) 步近似解出(构造在每簇内有零点的多项式)。图 5.5:\(n=14\),四簇(140、120 各一个,10 附近 10 个,其余在 \([0.95,1.05]\)),4 步后误差已很小,6 步后精确。

基于条件数的粗略界(\(\kappa(A)=\|A\|_2\|A^{-1}\|_2=\lambda_n/\lambda_1\);原书误印为 \(\lambda_1/\lambda_n\)):

\[\|x_k-x^*\|_A\le2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^k\|x_0-x^*\|_A. \tag{5.35}\]
(PDF 抽取文本作 \((\cdot)^{2k}\) 且无系数 2;标准结果为 \(2(\cdot)^k\),编写时宜采用标准形式并注明。)常严重高估误差,但只知道极端特征值时有用。与最速下降 (3.29) 形式相同,但依赖 \(\sqrt\kappa\) 而非 \(\kappa\)——这是 CG 比最速下降快得多的根本原因。例:\(\kappa=10^4\) 时最速下降每步缩减因子约 \(1-2\times10^{-4}\),CG 约 \(1-0.02\)。

预条件(Preconditioning,PDF p.118–119)

变量变换 \(\hat x=Cx\)(5.36),\(C\) 非奇异:

\[\hat\phi(\hat x)=\tfrac12\hat x^T(C^{-T}AC^{-1})\hat x-(C^{-T}b)^T\hat x, \tag{5.37}\]
即求解 \((C^{-T}AC^{-1})\hat x=C^{-T}b\),收敛速度取决于 \(C^{-T}AC^{-1}\) 的特征值。目标:选 \(C\) 使其条件数远小于 \(\kappa(A)\),或特征值成簇。无需显式变换,推导后算法只用到 \(M=C^TC\)(对称正定)。

Algorithm 5.3(预条件 CG):

给定 x_0,预条件子 M;r_0 ← A x_0 − b;解 M y_0 = r_0;p_0 ← −y_0;k ← 0
while r_k ≠ 0
  α_k ← r_kᵀy_k / p_kᵀA p_k               (5.38a)
  x_{k+1} ← x_k + α_k p_k                  (5.38b)
  r_{k+1} ← r_k + α_k A p_k                (5.38c)
  解 M y_{k+1} = r_{k+1}                    (5.38d)
  β_{k+1} ← r_{k+1}ᵀy_{k+1} / r_kᵀy_k      (5.38e)
  p_{k+1} ← −y_{k+1} + β_{k+1} p_k         (5.38f)
end

(PDF 文本初始化写作 \(p_0=-r_0\),应为 \(p_0=-y_0\)。)\(M=I\) 退化为标准 CG。残差正交性变为

\[r_i^TM^{-1}r_j=0,\quad i\neq j. \tag{5.39}\]
与无预条件相比的主要额外代价是每步解 \(My=r\)。

实用预条件子(Practical Preconditioners,PDF p.119–120)

  • 没有对所有矩阵都“最好”的预条件:需在 \(M\) 的有效性、构造/存储代价、解 \(My=r\) 的代价间权衡。
  • 特定类型矩阵(如偏微分方程离散化)有好策略:\(My=r\) 常是原系统的简化版(如更粗网格离散)。了解问题结构与来源是设计有效方法的关键。
  • 通用预条件子:对称逐次超松弛(SSOR)、不完全 Cholesky(incomplete Cholesky)、带状预条件。不完全 Cholesky 通常最有效:按 Cholesky 过程计算但只保留比真因子 \(L\) 更稀疏的近似 \(\tilde L\)(通常不比 \(A\) 的下三角更稠),\(A\approx\tilde L\tilde L^T\),取 \(C=\tilde L^T\),\(M=\tilde L\tilde L^T\),\(C^{-T}AC^{-1}=\tilde L^{-1}A\tilde L^{-T}\approx I\)。不显式算 \(M\),存 \(\tilde L\),解 \(My=r\) 用两次三角回代,代价与 \(Ap\) 相当。
  • 陷阱:结果可能不(够)正定,需增大对角元;稀疏性限制可能导致数值不稳定或中断,可允许更多填充,但更贵。

5.2 非线性共轭梯度法(Nonlinear Conjugate Gradient Methods,PDF p.140–151)

Fletcher–Reeves 方法(PDF p.140–141)

对 Algorithm 5.2 做两处改动:(1) 步长 \(\alpha_k\) 改为沿 \(p_k\) 的线搜索近似极小;(2) 残差 \(r\) 换成非线性目标的梯度。

Algorithm 5.4(FR-CG):

给定 x_0;计算 f_0, ∇f_0;p_0 ← −∇f_0;k ← 0
while ∇f_k ≠ 0
  计算 α_k,x_{k+1} = x_k + α_k p_k;计算 ∇f_{k+1}
  β^{FR}_{k+1} ← ∇f_{k+1}ᵀ∇f_{k+1} / ∇f_kᵀ∇f_k     (5.40a)
  p_{k+1} ← −∇f_{k+1} + β^{FR}_{k+1} p_k           (5.40b)
end
  • \(f\) 为强凸二次且精确线搜索时即线性 CG。每步只需函数值和梯度,无矩阵运算,只存几个向量,适合大规模。
  • 下降性需要线搜索条件:\(\nabla f_k^Tp_k=-\|\nabla f_k\|^2+\beta_k^{FR}\nabla f_k^Tp_{k-1}\)(5.41)。精确线搜索时 \(\nabla f_k^Tp_{k-1}=0\),必为下降方向;非精确时第二项可能占优使 \(p_k\) 成为上升方向。解决:步长满足强 Wolfe 条件
    \[f(x_k+\alpha_kp_k)\le f(x_k)+c_1\alpha_k\nabla f_k^Tp_k,\qquad |\nabla f(x_k+\alpha_kp_k)^Tp_k|\le c_2|\nabla f_k^Tp_k|, \tag{5.42}\]
    且 \(0<c_1<c_2<\frac12\)(注意比第 3 章的 \(c_2<1\) 更严)。由引理 5.6,这保证所有方向下降。

Polak–Ribière 方法(PDF p.141–142)

\[\beta_{k+1}^{PR}=\frac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{\|\nabla f_k\|^2}. \tag{5.43}\]
  • 强凸二次 + 精确线搜索时梯度正交 (5.15),\(\beta^{PR}=\beta^{FR}\);一般非线性 + 非精确线搜索时表现差异显著,数值经验表明 PR-CG 更稳健高效。
  • 意外事实:强 Wolfe 条件不能保证 PR 方向下降。定义
    \[\beta_{k+1}^+=\max\{\beta_{k+1}^{PR},0\}, \tag{5.44}\]
    得到 PR+ 算法,对强 Wolfe 稍作修改即可保证下降性。
  • Hestenes–Stiefel 公式:
    \[\beta_{k+1}^{HS}=\frac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{(\nabla f_{k+1}-\nabla f_k)^Tp_k}, \tag{5.45}\]
    理论与实践都与 PR 相似。推导:要求相邻方向关于线段 \([x_k,x_{k+1}]\) 上的平均 Hessian \(\bar G_k=\int_0^1\nabla^2f(x_k+\tau\alpha_kp_k)d\tau\) 共轭;由 \(\nabla f_{k+1}=\nabla f_k+\alpha_k\bar G_kp_k\),\(p_{k+1}^T\bar G_kp_k=0\) 即得 (5.45)。
  • 其他 \(\beta\) 选择都没有显著优于 PR。

二次终止与重启(Quadratic Termination and Restarts,PDF p.142–144)

  • 实现通常在线搜索中包含沿 \(p_k\) 的二次(或三次)插值,保证 \(f\) 为严格凸二次时步长精确,从而退化为线性 CG。
  • 重启:每 \(n\) 步令 \(\beta_k=0\)(走最速下降),周期性清除可能无益的旧信息。理论结果:\(n\) 步二次收敛
    \[\|x_{k+n}-x^*\|=O(\|x_k-x^*\|^2). \tag{5.46}\]
    直观:若 \(f\) 在解附近是强凸二次,迭代进入该区域后某次重启,之后就是线性 CG,\(n\) 步内终止;重启重要,因为线性 CG 的有限终止性要求 \(p_0\) 为负梯度。非精确二次时也有显著进展。
  • 但实践中意义有限:非线性 CG 只推荐用于大 \(n\),常在少于 \(n\) 步就得到近似解,重启可能从未发生。故常不重启,或用其他准则:最流行的是基于二次函数梯度正交性 (5.15),当相邻梯度远非正交时重启:
    \[\frac{|\nabla f_k^T\nabla f_{k-1}|}{\|\nabla f_k\|^2}\ge\nu,\quad \nu\approx0.1. \tag{5.47}\]
  • 其他:重启方向不用最速下降,如 Harwell VA14 用基于 \(\nabla f_{k+1}\)、\(p_k\) 及第三个历史方向的三项递推;CONMIN(第 9 章)更进一步。(5.44) 也可视为一种重启(\(\beta^{PR}<0\) 时回到最速下降),但因 \(\beta^{PR}\) 多数时候为正,这种重启很少发生。

数值表现(Numerical Performance,PDF p.144)

表 5.1:不重启的 FR、PR、PR+,强 Wolfe 参数 \(c_1=10^{-4}\)、\(c_2=0.1\),终止条件 \(\|\nabla f_k\|_\infty<10^{-5}(1+|f_k|)\) 或 10000 次迭代(记 *)。“mod”列为 PR+ 需要 (5.44) 调整的次数。

问题 n FR it/f-g PR it/f-g PR+ it/f-g mod
CALCVAR3 200 2808/5617 2631/5263 2631/5263 0
GENROS 500 * 1068/2151 1067/2149 1
XPOWSING 1000 533/1102 212/473 97/229 3
TRIDIA1 1000 264/531 262/527 262/527 0
MSQRT1 1000 422/849 113/231 113/231 0
XPOWELL 1000 568/1175 212/473 97/229 3
TRIGON 1000 231/467 40/92 40/92 0

FR 在 GENROS 上远离解时步长极短、改进极小。PR/PR+ 并非总比 FR 好,且多需一个向量存储,但作者推荐尽量用 PR-CG 或 PR+。

Fletcher–Reeves 方法的行为(PDF p.144–147)

假设水平集 \(\mathcal{L}\) 有界、\(f\) 二阶连续可微,由引理 3.1 存在强 Wolfe 步长。

引理 5.6:Algorithm 5.4 的步长满足强 Wolfe (5.42) 且 \(0<c_2<\frac12\),则所有 \(p_k\) 都是下降方向且

\[-\frac{1}{1-c_2}\le\frac{\nabla f_k^Tp_k}{\|\nabla f_k\|^2}\le\frac{2c_2-1}{1-c_2},\quad k=0,1,\dots \tag{5.48}\]
证明:\(t(\xi)=(2\xi-1)/(1-\xi)\) 在 \([0,\frac12]\) 单调增,\(t(0)=-1\),\(t(\frac12)=0\),故 \(-1<\frac{2c_2-1}{1-c_2}<0\)(5.49)。归纳:\(k=0\) 时中间项为 \(-1\) 满足。由 (5.40),
\[\frac{\nabla f_{k+1}^Tp_{k+1}}{\|\nabla f_{k+1}\|^2}=-1+\beta_{k+1}\frac{\nabla f_{k+1}^Tp_k}{\|\nabla f_{k+1}\|^2}=-1+\frac{\nabla f_{k+1}^Tp_k}{\|\nabla f_k\|^2}, \tag{5.50}\]
由 (5.42b) \(|\nabla f_{k+1}^Tp_k|\le-c_2\nabla f_k^Tp_k\),得 \(-1+c_2\frac{\nabla f_k^Tp_k}{\|\nabla f_k\|^2}\le\cdot\le-1-c_2\frac{\nabla f_k^Tp_k}{\|\nabla f_k\|^2}\),代入归纳假设左端得 \(-1-\frac{c_2}{1-c_2}\le\cdot\le-1+\frac{c_2}{1-c_2}\),即 (5.48) 对 \(k+1\) 成立。 只用到第二个强 Wolfe 条件;(5.48) 限制了 \(\|p_k\|\) 增长速度,在收敛分析中关键。

FR 的弱点:由 (5.48),\(\chi_1\frac{\|\nabla f_k\|}{\|p_k\|}\le\cos\theta_k\le\chi_2\frac{\|\nabla f_k\|}{\|p_k\|}\)(5.52)。若 \(p_k\) 很差(\(\cos\theta_k\approx0\)),则 \(\|\nabla f_k\|\ll\|p_k\|\);步长可能极小,\(x_{k+1}\approx x_k\),\(\nabla f_{k+1}\approx\nabla f_k\),于是 \(\beta_{k+1}^{FR}\approx1\)(5.53),\(p_{k+1}\approx p_k\)——新方向几乎不改进,将跟随一长串无效迭代。 PR 在此情形下:\(\nabla f_{k+1}\approx\nabla f_k\) 使 \(\beta^{PR}_{k+1}\approx0\),\(p_{k+1}\) 接近最速下降,\(\cos\theta_{k+1}\approx1\)——PR 遇坏方向后自动重启。PR+ 与 HS 同理。 实例(Gilbert–Nocedal):\(n=100\) 的问题上,\(\cos\theta_k\sim10^{-2}\) 持续数百步,步长 \(\sim10^{-2}\),FR 需数千步,PR 仅 37 步;FR 周期性沿最速下降重启会好很多。FR 不应在无重启策略下实现。

全局收敛(Global Convergence,PDF p.147–151)

与线性 CG 不同,非线性 CG 的收敛性质令人惊讶、有时古怪,理论仍不完整。

假设 5.1:(i) 水平集 \(\mathcal{L}=\{x:f(x)\le f(x_0)\}\) 有界;(ii) 在 \(\mathcal{L}\) 的某邻域 \(\mathcal{N}\) 内梯度 Lipschitz 连续(5.54)。由此存在 \(\bar\gamma\) 使 \(\|\nabla f(x)\|\le\bar\gamma\),\(x\in\mathcal{L}\)(5.55)。

定理 5.7(Zoutendijk 重述):在假设 5.1 下,下降方向 + Wolfe 步长,则 \(\sum_{k\ge1}\cos^2\theta_k\|\nabla f_k\|^2<\infty\)(5.56)。

  • 周期重启的版本:重启步 \(\cos\theta=1\),\(\sum_{k=k_1,k_2,\dots}\|\nabla f_k\|^2<\infty\)(5.57);若两次重启间隔不超过 \(\bar n\),重启步无穷多,得 \(\liminf\|\nabla f_k\|=0\)(5.58),适用于本章所有算法的重启版。
  • 更有意义的是不重启版本(大规模问题 \(n\ge1000\) 通常在远少于 \(n\) 步内收敛,重启从不发生)。

定理 5.8(FR 全局收敛):假设 5.1 成立,Algorithm 5.4 步长满足强 Wolfe 且 \(0<c_1<c_2<\frac12\),则

\[\liminf_{k\to\infty}\|\nabla f_k\|=0. \tag{5.59}\]
证明(反证,设 \(\|\nabla f_k\|\ge\gamma>0\)(5.60)):

  1. 由引理 5.6,\(\cos\theta_k\ge\frac{1-2c_2}{1-c_2}\frac{\|\nabla f_k\|}{\|p_k\|}\)(5.61,PDF 文本系数印作 \(\frac{1}{1-c_2}\),按 (5.48) 右端应为 \(\frac{1-2c_2}{1-c_2}\),均为正常数,不影响结论),代入 Zoutendijk 得 \(\sum\frac{\|\nabla f_k\|^4}{\|p_k\|^2}<\infty\)(5.62)。
  2. 由 (5.42b) 与引理 5.6,\(|\nabla f_k^Tp_{k-1}|\le-c_2\nabla f_{k-1}^Tp_{k-1}\le\frac{c_2}{1-c_2}\|\nabla f_{k-1}\|^2\)(5.63),故 \(\|p_k\|^2\le\|\nabla f_k\|^2+2\beta_k^{FR}|\nabla f_k^Tp_{k-1}|+(\beta_k^{FR})^2\|p_{k-1}\|^2\le c_3\|\nabla f_k\|^2+(\beta_k^{FR})^2\|p_{k-1}\|^2\),\(c_3=\frac{1+c_2}{1-c_2}\)。递推并利用 \((\beta_k^{FR})^2\cdots(\beta_{k-i}^{FR})^2=\|\nabla f_k\|^4/\|\nabla f_{k-i-1}\|^4\) 和 \(p_1=-\nabla f_1\):
    \[\|p_k\|^2\le c_3\|\nabla f_k\|^4\sum_{j=1}^k\|\nabla f_j\|^{-2}. \tag{5.64}\]
  3. 由 (5.55)(5.60),\(\|p_k\|^2\le\frac{c_3\bar\gamma^4}{\gamma^2}k\)(5.65),故 \(\sum\frac1{\|p_k\|^2}\ge\gamma_4\sum\frac1k=\infty\)(5.66)。
  4. 而由 (5.62) 与 (5.60),\(\sum\frac1{\|p_k\|^2}<\infty\)(5.67),矛盾。

此结果适用于 FR 的实用实现与一般非线性函数,比只对凸函数成立的结果更令人满意。

一般地,若存在 \(c_4,c_5>0\) 使 \(\cos\theta_k\ge c_4\frac{\|\nabla f_k\|}{\|p_k\|}\) 且 \(\frac{\|\nabla f_k\|}{\|p_k\|}\ge c_5\),则由定理 5.7 得 \(\lim\|\nabla f_k\|=0\)。对 PR 方法,可在 \(f\) 强凸 + 精确线搜索下证明。

对一般非凸函数,无法对 PR 证明类似定理 5.8 的结果(尽管实践中 PR 更好)。

定理 5.9(Powell 反例):PR 方法 + 理想线搜索(取 \(t(\alpha)=f(x_k+\alpha p_k)\) 的第一个驻点),存在二阶连续可微 \(f:\mathbb{R}^3\to\mathbb{R}\) 及初始点,使 \(\|\nabla f_k\|\) 远离零——PR 可能无限循环而不接近解。证明复杂,只证存在性不显式构造。所假设的步长(第一个驻点)可被任何实用线搜索接受。证明需要相邻搜索方向几乎互为相反,理想线搜索下只有 \(\beta_k<0\) 才可能,这启发了 \(\beta_k^+=\max\{\beta_k^{PR},0\}\)(5.68),即 PR+。配合修改的 Wolfe 线搜索保证下降,PR+ 对一般函数全局收敛。

注释与参考(PDF p.151–152)

  • CG 1950 年代由 Hestenes–Stiefel 提出作为对称正定系统精确解法;多年后才被视为能在远少于 \(n\) 步给出好近似的迭代法(稀疏线性代数最重要进展之一)。本章线性 CG 讲法沿用 Luenberger。
  • FR 非线性 CG 提出于线性 CG 失宠之后、被重新发现为迭代法之前。PR 见 Polak–Ribière;非凸失败反例见 Powell;重启见 Powell。
  • Powell 进一步分析 FR + 精确线搜索的低效:若迭代进入 \(f=\frac12x^Tx\) 的二维区域,梯度与搜索方向夹角保持不变,若接近 90° 则收敛极慢,可比最速下降还慢;PR 则遇小步后趋于最速下降,避免一串小步。全局收敛另见 Al-Baali、Gilbert–Nocedal。
  • 收敛速度(多假设精确线搜索):Crowder–Wolfe 证明线性收敛且构造例子说明不能 Q-超线性;Powell:进入二次区域则要么有限终止要么线性;Cohen、Burmeister:一般函数 \(n\) 步二次收敛;Ritter:实际为 \(n\) 步超二次 \(o(\|x_k-x^*\|^2)\);方向一致线性无关时更快(难验证)。
  • Nemirovsky–Yudin:强凸问题上 FR、PR 达不到最优复杂度界,甚至可能比最速下降慢;Nesterov 提出达到最优界的算法(与 PARTAN 平行切线法相关),作者当时认为不太可能实用(注:这就是后来广泛使用的 Nesterov 加速梯度法,此为 1999 年观点)。
  • PR 保证全局收敛的特殊线搜索有缺点;表 5.1 来自 Gilbert–Nocedal,该文给出保证 PR+ 下降的线搜索并证明全局收敛。

习题概览(PDF p.152–153)

  • 5.1 实现 CG 解 Hilbert 矩阵 \(A_{ij}=1/(i+j-1)\) 方程组,\(b=(1,\dots,1)^T\),\(x_0=0\),\(n=5,8,12,20\),报告残差降到 \(10^{-6}\) 的迭代次数(体会病态矩阵)。
  • 5.2 共轭向量线性无关;5.3 验证 (5.6);5.4 \(h(\sigma)\) 严格凸;5.5 直接验证 \(k=1\) 时 (5.16)(5.17);5.6 (5.23d) ⇔ (5.13d)。
  • 5.7 \([I+P_k(A)A]^TA[I+P_k(A)A]\) 的特征对为 \(\lambda_i[1+\lambda_iP_k(\lambda_i)]^2\)、\(v_i\)。
  • 5.8 构造不同特征值分布的矩阵检验定理 5.5;5.9 推导预条件 CG;5.10 验证 (5.39)。
  • 5.11 二次函数 + 精确线搜索时 PR、HS 退化为 FR;5.12 引理 5.6 对任何 \(|\beta_k|\le\beta_k^{FR}\) 成立。

本章要点

  1. 共轭方向:关于 \(A\) 共轭的方向上逐次精确极小,\(n\) 步终止;扩展子空间极小化(残差与之前方向正交)。
  2. CG 只用前一方向即自动共轭;残差两两正交;方向与残差张成 Krylov 子空间。
  3. 实用 CG(Algorithm 5.2)每步一次矩阵–向量乘;适合大规模稀疏问题。
  4. 收敛速度由特征值分布决定:\(r\) 个不同特征值 \(r\) 步终止;成簇则快;条件数界 \(\left(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\right)\) 优于最速下降的 \(\frac{\kappa-1}{\kappa+1}\)。
  5. 预条件:改善 \(C^{-T}AC^{-1}\) 的谱,不完全 Cholesky 最常用。
  6. 非线性 CG:FR、PR、PR+、HS;强 Wolfe 且 \(c_2<\frac12\) 保证 FR 下降与 \(\liminf\|\nabla f_k\|=0\);FR 遇坏方向陷入低效,PR 自动重启;PR 非凸可能不收敛(Powell 反例),PR+ 全局收敛;推荐 PR/PR+。

与量化交易的关联

  • 大规模线性方程组:风险模型中求解 \(\Sigma w=\mu\)(最大夏普/最小方差组合的无约束解)、广义最小二乘、岭回归正规方程。资产数达数千、使用“因子 + 特质”结构协方差 \(\Sigma=BFB^T+D\) 时,矩阵–向量乘只需 \(O(nK)\),可用 CG 免去显式构造 \(n\times n\) 矩阵;更简单的做法是用 Woodbury 公式,但 CG 在加入其他结构(如交易成本二次项)时更通用。
  • 特征值聚类的直接应用:因子协方差矩阵恰有“少数大特征值(\(K\) 个因子)+ 大量相近的小特征值(特质方差)”的结构,正是图 5.3 的情形——CG 约 \(K+1\) 步就能得到好解。用对角特质方差 \(D\) 做预条件(Jacobi 预条件)后,\(D^{-1/2}\Sigma D^{-1/2}=I+\) 低秩,特征值严格成两簇,CG 收敛极快。这是在大规模组合优化中用迭代法的有力论据。
  • 非线性 CG 用于大规模校准/机器学习:参数众多的模型(如大规模因子模型的联合估计)在内存受限时可用 PR+ CG;但现代实践中更多用 L-BFGS(第 9 章)。SciPy 的 method='CG' 就是 PR 变体。
  • Nesterov 加速:本章注释提到的 Nesterov 方法后来成为大规模凸优化(如组合优化的一阶方法、近端梯度/FISTA 求解 Lasso 型稀疏组合)的主力,可作为教材扩展。
  • 小规模问题(几十到几百资产)直接用 Cholesky 分解更稳健,本章明确指出 CG 只推荐用于大规模问题。

推荐习题

  • 5.1(Hilbert 矩阵上的 CG,体会病态与舍入误差)。
  • 5.8(构造特征值分布,验证聚类效应——可直接用“因子 + 特质”结构协方差做实验)。
  • 5.9、5.10(推导预条件 CG)。
  • 5.11(理解 FR/PR/HS 的联系)。
  • 5.12(引理 5.6 的推广,理解 \(\beta\) 的取值范围与下降性)。

第 6 章 实用牛顿法(Practical Newton Methods,PDF p.154–183)

6.0 引言(PDF p.155–156)

  • 纯牛顿法(单位步)接近 \(x^*\) 时收敛快,但从远处起点可能不收敛,在非凸区域行为不稳定。目标:在所有情形下稳健且高效。
  • 牛顿步解对称线性方程组
    \[\nabla^2f(x_k)p_k^N=-\nabla f(x_k). \tag{6.1}\]
    Hessian 正定时 \(p_k^N\) 是下降方向;不正定或接近奇异时,可能是上升方向或过长。
  • 两种保证步长质量的策略:(1) 用 CG 解 (6.1),遇负曲率即终止——Newton–CG(有线搜索与信赖域两种实现);(2) 在解 (6.1) 之前或过程中修改 Hessian 使其充分正定——修正牛顿法(modified Newton)。
  • 控制计算量:Newton–CG 在得到精确解前就停止 CG,即非精确牛顿(inexact Newton);直接法则利用 Hessian 稀疏结构做稀疏高斯消元。
  • Hessian 计算通常是主要工作;无解析式时用自动微分或有限差分(第 7 章)。综合修正、稀疏性与微分技术,本章的牛顿法是求解中小及大规模无约束问题最可靠强大的方法之一。

6.1 非精确牛顿步(Inexact Newton Steps,PDF p.156–158)

  • 大规模时分解代价高;且远离解时二次模型本就不准,精确求解 (6.1) 不值得。用迭代法求近似解。
  • 终止规则基于残差
    \[r_k=\nabla^2f(x_k)p_k+\nabla f(x_k), \tag{6.2}\]
    因 \(r\) 对目标的数乘缩放不具不变性,取相对于右端的大小:
    \[\|r_k\|\le\eta_k\|\nabla f(x_k)\|, \tag{6.3}\]
    \(\{\eta_k\}\)(\(0<\eta_k<1\))称为强制序列(forcing sequence)。

定理 6.1:设 \(\nabla f\) 在极小点 \(x^*\) 邻域内连续可微,\(\nabla^2f(x^*)\) 正定。迭代 \(x_{k+1}=x_k+p_k\),\(p_k\) 满足 (6.3),且 \(\eta_k\le\eta\in[0,1)\)。若 \(x_0\) 足够接近 \(x^*\),则 \(\{x_k\}\) 线性收敛:\(\|x_{k+1}-x^*\|\le c\|x_k-x^*\|\),\(0<c<1\)。

  • 条件不苛刻:若允许 \(\eta_k\ge1\),\(p_k=0\) 总满足 (6.3),迭代停滞。
  • 非正式推导:\(\|\nabla^2f(x_k)^{-1}\|\le L\),故 \(\|p_k\|\le L(\|\nabla f_k\|+\|r_k\|)\le2L\|\nabla f_k\|\);由 Taylor,
    \[\nabla f(x_{k+1})=\nabla f_k+\nabla^2f_kp_k+O(\|p_k\|^2)=r_k+O(\|\nabla f_k\|^2), \tag{6.4}\]
    \[\|\nabla f(x_{k+1})\|\le\eta_k\|\nabla f_k\|+O(\|\nabla f_k\|^2). \tag{6.5}\]
    故梯度每步约缩减 \(\eta_k\) 倍,\(\limsup\frac{\|\nabla f_{k+1}\|}{\|\nabla f_k\|}\le\eta<1\)。若 \(r_k=o(\|\nabla f_k\|)\) 则超线性;若 \(r_k=O(\|\nabla f_k\|^2)\) 则 \(\limsup\frac{\|\nabla f_{k+1}\|}{\|\nabla f_k\|^2}=c\),恢复二次收敛。迭代点与梯度以相同速度收敛。

定理 6.2:在定理 6.1 条件下且 \(x_k\to x^*\),若 \(\eta_k\to0\) 则超线性收敛;若 \(\eta_k=O(\|\nabla f(x_k)\|)\) 则二次收敛(习题 6.6)。

  • 例:\(\eta_k=\min(0.5,\sqrt{\|\nabla f_k\|})\) 得超线性;\(\eta_k=\min(0.5,\|\nabla f_k\|)\) 得二次。
  • 这些结果(Dembo–Eisenstat–Steihaug)都是局部的,假设迭代已进入解附近且取单位步(全局化策略不妨碍快速收敛)。以下各节说明非精确牛顿可嵌入实用的线搜索与信赖域实现。

6.2 线搜索牛顿法(Line Search Newton Methods,PDF p.159–162)

\(x_{k+1}=x_k+\alpha_kp_k\),\(p_k\) 为牛顿方向或其近似;\(\alpha_k\) 满足 Wolfe/Goldstein 或 Armijo 回溯,总是先试单位步。

线搜索 Newton–CG(Line Search Newton–CG,PDF p.159–161)

又称截断牛顿法(truncated Newton):用 CG 解 (6.1) 并尝试满足 (6.3)。由于 CG 为正定系统设计,而远离解时 Hessian 可能有负特征值,一旦生成负曲率方向就终止 CG——保证 \(p_k\) 下降且保持牛顿法快速收敛。内层 CG(\(A=\nabla^2f_k\),\(b=-\nabla f_k\),上标 \((i)\) 表示 CG 内部迭代)三个要求:

  • (a) CG 起点 \(x^{(0)}=0\);
  • (b) 负曲率检验:若方向满足 \((p^{(i)})^TAp^{(i)}\le0\)(6.6),若是第一次 CG 迭代(\(i=0\)),完成该步算出 \(x^{(1)}\) 后停止;若 \(i>0\),立即停止并返回 \(x^{(i)}\);
  • (c) 牛顿步 \(p_k\) 取最终 CG 迭代点 \(x^{(f)}\)。

说明:也可取 \(x^{(0)}\) 为上一个牛顿–CG 步,两种选择在线搜索中表现相当。若首步就遇负曲率,\(p_k\) 即最速下降方向 \(-\nabla f_k\)(这正是取 \(x^{(0)}=0\) 的原因);若 CG 执行多于一步,\(p_k\) 也是下降方向(习题 6.2)。可在 CG 中用预条件。

Hessian-free:Newton–CG 不需显式 Hessian,只需 Hessian–向量积 \(\nabla^2f(x_k)p\)。有限差分:

\[\nabla^2f(x_k)p\approx\frac{\nabla f(x_k+hp)-\nabla f(x_k)}{h}, \tag{6.7}\]
精度 \(O(h)\),每次 CG 迭代多一次梯度计算;或用自动微分精确计算。这类方法称 Hessian-free 牛顿法。

Algorithm 6.1(线搜索 Newton–CG):

给定 x_0
for k = 0,1,2,...
  对 ∇²f(x_k) p = −∇f_k 从 x^(0)=0 开始做 CG;
  当 ‖r_k‖ ≤ min(0.5, √‖∇f_k‖)‖∇f_k‖ 或遇负曲率(按 (b))时终止
  x_{k+1} = x_k + α_k p_k,α_k 满足 Wolfe、Goldstein 或 Armijo 回溯
end

(正文中强制序列印作 \(\min(0.5,\sqrt{\|\nabla^2f_k\|})\),算法中为 \(\sqrt{\|\nabla f_k\|}\),后者正确。)

弱点(尤其无预条件时):Hessian 接近奇异时 Newton–CG 方向可能过长,线搜索需多次函数求值,且下降很小。规范化牛顿步的好规则难定(可能破坏良好尺度时的快速收敛);在 (6.6) 中引入阈值更好但阈值也难选。6.4 节的信赖域 Newton–CG 处理得更好,作者略偏好后者。

修正牛顿法(Modified Newton's Method,PDF p.161–162)

希望用直接法(如高斯消元)解 (6.1)。Hessian 不正定或接近奇异时,在求解前或过程中修改它(加正对角阵或满矩阵),得到正定近似。

Algorithm 6.2(带修正的线搜索牛顿法):

给定 x_0
for k = 0,1,2,...
  分解 B_k = ∇²f(x_k) + E_k;∇²f(x_k) 充分正定时 E_k = 0,否则选 E_k 使 B_k 充分正定
  解 B_k p_k = −∇f(x_k)
  x_{k+1} = x_k + α_k p_k,α_k 满足 Wolfe、Goldstein 或 Armijo 回溯
end

\(E_k\) 的选择至关重要;有些方法不显式算 \(E_k\),而在标准分解过程中“边分解边修改”。

有界修正分解性质(bounded modified factorization property):只要 Hessian 序列有界,\(B_k\) 条件数就一致有界:

\[\mathrm{cond}(B_k)=\|B_k\|\|B_k^{-1}\|\le C. \tag{6.8}\]

定理 6.3:\(f\) 在开集 \(\mathcal{D}\) 上二阶连续可微,水平集 \(\mathcal{L}=\{x\in\mathcal{D}:f(x)\le f(x_0)\}\) 紧,有界修正分解性质成立,则 \(\lim\nabla f(x_k)=0\)。 证明:线搜索保证迭代留在水平集;Hessian 连续 + 紧 ⇒ Hessian 有界 ⇒ (6.8);由 Zoutendijk 定理及 (3.16)(\(\cos\theta_k\ge1/C\))得证。

收敛速度:若 \(x_k\to x^*\) 且 \(\nabla^2f(x^*)\) 充分正定使修正最终为零,由定理 3.5 最终 \(\alpha_k=1\),退化为纯牛顿,二次收敛。若 \(\nabla^2f^*\) 接近奇异,修正可能不消失,只能线性收敛。

6.3 Hessian 修正(Hessian Modifications,PDF p.162–174)

目标:\(B_k=\nabla^2f(x_k)+E_k\) 充分正定且条件良好(使定理 6.3 成立);修正尽量小以保留二阶信息;分解代价适中。先讲基于特征分解的“理想”策略。

特征值修正(Eigenvalue Modification,PDF p.163–164)

例:\(\nabla f(x_k)=(1,-3,2)\),\(\nabla^2f(x_k)=\mathrm{diag}(10,3,-1)\)(不定),\(Q=I\),

\[\nabla^2f(x_k)=Q\Lambda Q^T=\sum_{i=1}^n\lambda_iq_iq_i^T. \tag{6.9}\]
纯牛顿步 \(p_k^N=(-0.1,1,2)\),\(\nabla f_k^Tp_k^N>0\),不是下降方向。

  • 把负特征值换成小正数 \(\delta=\sqrt u\)(\(u=10^{-16}\) 为机器精度):
    \[B_k=\sum_{i=1}^2\lambda_iq_iq_i^T+\delta q_3q_3^T=\mathrm{diag}(10,3,10^{-8}), \tag{6.10}\]
    保留了 \(q_1,q_2\) 方向的曲率,但
    \[p_k=-B_k^{-1}\nabla f_k=-\sum_{i=1}^2\frac1{\lambda_i}q_i(q_i^T\nabla f_k)-\frac1\delta q_3(q_3^T\nabla f(x_k))\approx-(2\times10^8)q_3, \tag{6.11}\]
    几乎平行于 \(q_3\) 且极长。虽然 \(f\) 沿此方向下降,但极端长度违背牛顿法依赖局部二次模型的精神,有效性存疑。
  • 其他选择:翻转负特征值符号(\(\delta=1\));令 (6.11) 末项为零(不含负曲率分量);按步长不过长自适应选 \(\delta\)(有信赖域味道)。哪种最理想尚无共识。
  • 最小 Frobenius 范数修正:对称 \(A=Q\Lambda Q^T\),使 \(\lambda_{\min}(A+\Delta A)\ge\delta\) 的最小 Frobenius 范数修正为
    \[\Delta A=Q\,\mathrm{diag}(\tau_i)Q^T,\quad \tau_i=\begin{cases}0,&\lambda_i\ge\delta,\\\delta-\lambda_i,&\lambda_i<\delta.\end{cases} \tag{6.12}\]
    (\(\|A\|_F^2=\sum a_{ij}^2\))一般非对角,\(A+\Delta A=Q(\Lambda+\mathrm{diag}(\tau_i))Q^T\)。(6.10) 就是这种最优修正。
  • 最小欧氏(2-)范数修正:
    \[\Delta A=\tau I,\quad \tau=\max(0,\delta-\lambda_{\min}(A)), \tag{6.13}\]
    修正后为 \(A+\tau I\)(6.14),与(无尺度)信赖域中的矩阵形式相同,所有特征值都被平移到 \(\ge\delta\)。
  • 实用软件不用特征分解(太贵),而用高斯消元间接选择修正,数值经验表明常(但非总)能得到好方向。

加单位阵倍数(Adding a Multiple of the Identity,PDF p.164–165)

找 \(\tau>0\) 使 \(\nabla^2f(x_k)+\tau I\) 充分正定;由 (6.13) 需知最小特征值,通常无好估计。利用最大绝对特征值不超过 \(\|A\|_F\):

Algorithm 6.3(加单位阵倍数的 Cholesky):

β ← ‖A‖_F
if min_i a_ii > 0:  τ_0 ← 0   else  τ_0 ← β/2
for k = 0,1,2,...
  尝试 Cholesky 分解 L Lᵀ = A + τ_k I     (原文写“incomplete Cholesky”,此处应理解为对 A+τI 做 Cholesky)
  if 成功: 返回 L
  else: τ_{k+1} ← max(2τ_k, β/2)
end

简单,可能比后面的修正分解更可取;缺点:\(\tau\) 可能不必要地大,使方向过度偏向最速下降;每个 \(\tau_k\) 都需重新做数值分解(符号分解只做一次),试多个值时代价高。

修正 Cholesky 分解(Modified Cholesky Factorization,PDF p.165–170)

思路:对 \(\nabla^2f(x_k)\) 做 Cholesky,过程中必要时增大对角元。两个目标:保证修正因子存在且相对 Hessian 范数有界;Hessian 充分正定时不修改。

\(LDL^T\) 形式:对称正定 \(A\) 可写成

\[A=LDL^T, \tag{6.15}\]
\(L\) 单位下三角,\(D\) 对角且元素为正。

例 6.1(n=3):比较各列元素:\(d_1=a_{11}\),\(l_{21}=a_{21}/d_1\),\(l_{31}=a_{31}/d_1\);\(d_2=a_{22}-d_1l_{21}^2\),\(l_{32}=(a_{32}-d_1l_{31}l_{21})/d_2\);\(d_3=a_{33}-d_1l_{31}^2-d_2l_{32}^2\)。

Algorithm 6.4(Cholesky,\(LDL^T\) 形式):

for j = 1..n
  c_jj ← a_jj − Σ_{s=1}^{j−1} d_s l_js²;  d_j ← c_jj
  for i = j+1..n
    c_ij ← a_ij − Σ_{s=1}^{j−1} d_s l_is l_js;  l_ij ← c_ij / d_j
end

\(A\) 正定时所有 \(d_j>0\)。与标准形式 \(A=MM^T\)(6.16)的关系:\(M=LD^{1/2}\)。\(A\) 不定时 \(LDL^T\) 可能不存在;即使存在也数值不稳定(元素可任意大),因此“先分解后修改对角元”可能失败或得到与 \(A\) 差别极大的矩阵。

改为在分解过程中修改:选参数 \(\delta,\beta>0\),要求计算第 \(j\) 列时

\[d_j\ge\delta,\quad |m_{ij}|\le\beta,\ i=j+1,\dots,n, \tag{6.17}\]
(\(m_{ij}=l_{ij}\sqrt{d_j}\) 为 \(M\) 的元素)只需把 \(d_j\) 的计算改为
\[d_j=\max\left(|c_{jj}|,\left(\frac{\theta_j}{\beta}\right)^2,\delta\right),\quad \theta_j=\max_{j<i\le n}|c_{ij}|. \tag{6.18}\]
验证:\(c_{ij}=l_{ij}d_j\),\(|m_{ij}|=|l_{ij}\sqrt{d_j}|=\frac{|c_{ij}|}{\sqrt{d_j}}\le\frac{|c_{ij}|\beta}{\theta_j}\le\beta\)。\(\theta_j\) 可先于 \(d_j\) 算出(\(c_{ij}\) 不依赖 \(d_j\)),这正是引入 \(c_{ij}\) 的原因;计算应重排为先算 \(c_{ij}\) 再算 \(d_j\)。为减小修正量,引入对称行列交换,使第 \(j\) 步的行列具有最大对角元。

Algorithm 6.5(修正 Cholesky,Gill–Murray–Wright):

给定 δ>0, β>0
for k=1..n: c_kk ← a_kk                      (初始化对角元)
找 q 使 |c_qq| ≥ |c_ii|, i=j..n;交换第 j 与 q 行列
for j = 1..n                                  (计算 L 的第 j 列)
  for s = 1..j−1: l_js ← c_js / d_s
  for i = j+1..n: c_ij ← a_ij − Σ_{s=1}^{j−1} l_js c_is
  θ_j ← 0;if j < n: θ_j ← max_{j<i≤n} |c_ij|
  d_j ← max{|c_jj|, (θ_j/β)², δ}
  if j < n: for i = j+1..n: c_ii ← c_ii − c_ij² / d_j
end

(按文字说明,选主元交换应在每个 \(j\) 步进行。)约 \(n^3/6\) 次运算,与标准 Cholesky 相当,但行列交换的数据移动在大问题上可能代价可观;无额外存储(\(L,D,c_{ij}\) 可覆盖 \(A\))。

结果:设 \(P\) 为置换矩阵,

\[PAP^T+E=LDL^T=MM^T, \tag{6.19}\]
\(E\) 为非负对角阵,\(A\) 充分正定时为零;\(e_j=d_j-c_{jj}\),增大 \(c_{jj}\) 等价于增大原数据 \(a_{jj}\)。 参数:\(\delta=u\max(\gamma(A)+\xi(A),1)\),\(\gamma=\max_i|a_{ii}|\),\(\xi=\max_{i\neq j}|a_{ij}|\);Gill–Murray–Wright 建议
\[\beta=\max\left(\gamma(A),\frac{\xi(A)}{\sqrt{n^2-1}},u\right)^{1/2},\]
旨在最小化 \(\|E\|_\infty\)。

例 6.2:\(A=\begin{bmatrix}4&2&1\\2&6&3\\1&3&-0.004\end{bmatrix}\),特征值约 \(-1.25,2.87,8.38\)。Algorithm 6.5 给出 \(M=\begin{bmatrix}0.8165&1.8257&0\\2.4495&0&0\\1.2247&-1.2\times10^{-16}&1.2264\end{bmatrix}\)(含置换),\(E=\mathrm{diag}(0,0,3.008)\);修正后 \(A'=\begin{bmatrix}4&2&1\\2&6&3\\1&3&3.004\end{bmatrix}\),特征值 \(1.13,3.00,8.87\),条件数 7.8,相当温和。 Moré–Sorensen 证明:对精确 Hessian 用 Algorithm 6.5 得到的 \(B_k\) 条件数有界,(6.8) 成立。

Gershgorin 修正(Gershgorin Modification,PDF p.170)

先对 \(A\) 用 Algorithm 6.5 得 (6.19);若 \(E=0\) 结束。否则计算 \(\lambda_{\min}(A)\) 的两个上界:\(b_1\) 由 Gershgorin 圆盘定理(保证 \(A+b_1I\) 严格对角占优),\(b_2=\max_ie_{ii}\);取 \(\mu=\min(b_1,b_2)\),对 \(A+\mu I\) 做 Cholesky 作为修正因子。\(b_2\) 对较松的 \(b_1\) 起控制作用(Schnabel–Eskow)。是否优于单用 Algorithm 6.5 不明。两者都不修改“充分正定”的矩阵,但“充分”难以用 \(\lambda_{\min}\) 量化,可能修改 \(\lambda_{\min}>\delta\) 的矩阵。

修正对称不定分解(Modified Symmetric Indefinite Factorization,PDF p.171–174)

任何对称 \(A\) 可写成

\[PAP^T=LBL^T, \tag{6.20}\]
\(L\) 单位下三角,\(B\) 为 \(1\times1\) 与 \(2\times2\) 块的块对角,\(P\) 置换。允许 \(2\times2\) 块保证分解总存在且数值稳定(对不定矩阵算 \(LDL^T\) 不明智,因元素可能放大舍入误差)。

例 6.3:\(A=\begin{bmatrix}0&1&2&3\\1&2&2&2\\2&2&3&3\\3&2&3&4\end{bmatrix}\),\(P=[e_1,e_4,e_3,e_2]\), \(L=\begin{bmatrix}1&0&0&0\\0&1&0&0\\\frac19&\frac23&1&0\\\frac29&\frac13&0&1\end{bmatrix}\),\(B=\begin{bmatrix}0&3&0&0\\3&4&0&0\\0&0&\frac79&\frac59\\0&0&\frac59&\frac{10}9\end{bmatrix}\)(6.21),两个对角块都是 \(2\times2\)。

  • 惯性(inertia,正/零/负特征值个数):\(B\) 与 \(A\) 惯性相同;\(2\times2\) 块总构造成一正一负特征值,故 \(A\) 的正特征值数 = 正 \(1\times1\) 块数 + \(2\times2\) 块数。
  • 过程:选非奇异主元块 \(E\)(单个对角元或两个对角元加对应非对角元),置换使之成为前主子阵:\(\Pi A\Pi^T=\begin{bmatrix}E&C^T\\C&H\end{bmatrix}\)(6.22),块分解
    \[\Pi A\Pi^T=\begin{bmatrix}I&0\\CE^{-1}&I\end{bmatrix}\begin{bmatrix}E&0\\0&H-CE^{-1}C^T\end{bmatrix}\begin{bmatrix}I&E^{-1}C^T\\0&I\end{bmatrix},\]
    对 Schur 补 \(H-CE^{-1}C^T\) 递归。
  • 约 \(n^3/3\) 次浮点运算(文中称与正定 Cholesky 相同),外加选主元和置换的代价(可能可观)。理想的选主元策略:便宜、每步剩余矩阵元素增长有限、填充不过多。
  • Bunch–Parlett:搜索整个工作矩阵,找最大对角元 \(\xi_{dia}\) 与最大非对角元 \(\xi_{off}\);以对角元作 \(1\times1\) 主元时增长受 \(\xi_{dia}/\xi_{off}\) 控制,可接受则用之,否则取含 \(\xi_{off}\) 的 \(2\times2\) 块 \(\begin{bmatrix}a_{ii}&a_{ij}\\a_{ij}&a_{jj}\end{bmatrix}\)。数值稳定,\(L\) 最大元不超过 2.781;缺点是总共 \(O(n^3)\) 次比较,总时间可能约为 Algorithm 6.5 的两倍。
  • Bunch–Kaufman:每步至多搜索两列,总 \(O(n^2)\) 次比较;但 \(L\) 元素可能任意大,不适合修正 Cholesky 策略。
  • 有界 Bunch–Kaufman:折中——监控 \(L\) 元素大小,增长温和时用便宜的 BK 选择,过大时进一步搜索;通常代价接近 BK,最坏接近 BP。
  • 大型稀疏矩阵还须考虑主元选择对 \(L\) 稀疏性的影响(Duff 等、Duff–Reid、Fourer–Mehrotra)。
  • 修正:算出 (6.20) 后,对块对角 \(B\) 做谱分解 \(B=Q\Lambda Q^T\)(很便宜,逐块做,习题 6.7),构造
    \[F=Q\,\mathrm{diag}(\tau_i)Q^T,\quad \tau_i=\begin{cases}0,&\lambda_i\ge\delta,\\\delta-\lambda_i,&\lambda_i<\delta,\end{cases} \tag{6.23}\]
    即使 \(B+F\) 所有特征值不小于 \(\delta\) 的最小 Frobenius 修正。于是 \(P(A+E)P^T=L(B+F)L^T\),\(E=P^TLFL^TP\)(一般非对角)。与修正 Cholesky 不同,此法改变整个 \(A\) 而非仅对角;目标是 \(\lambda_{\min}(A)<\delta\) 时 \(\lambda_{\min}(A+E)\approx\delta\),但不一定总能接近(Cheng–Higham)。

6.4 信赖域牛顿法(Trust-Region Newton Methods,PDF p.174–182)

信赖域不要求模型 Hessian 正定,可直接用 \(B_k=\nabla^2f(x_k)\):

\[\min_p m_k(p)=f_k+\nabla f_k^Tp+\tfrac12p^TB_kp\quad\text{s.t. }\|p\|\le\Delta_k. \tag{6.24}\]
关注大规模实现,讨论四种技术:dogleg、二维子空间、迭代精确解、CG(并讨论信赖域范数的选择——相当于预条件)。

Newton–Dogleg 与子空间极小化(PDF p.174–175)

  • Hessian 不定时 dogleg 不直接适用;用 6.3 节的修正 Hessian 作 \(B_k\),保证凸二次,然后沿 dogleg 路径
    \[\tilde p(\tau)=\begin{cases}\tau p^U,&0\le\tau\le1,\\p^U+(\tau-1)(p^B-p^U),&1\le\tau\le2,\end{cases} \tag{6.25}\]
    极小化。二维子空间:
    \[\min_p m_k(p)\quad\text{s.t. }\|p\|\le\Delta_k,\ p\in\mathrm{span}\{\nabla f_k,p^B\}. \tag{6.26}\]
  • 优点:线性代数可全用直接法;全局收敛;希望修正分解在 Hessian 充分正定时不修改,保持局部快速收敛。
  • 不足:修正分解以“近乎随机”的方式扰动 Hessian,偏重某些方向,信赖域的好处可能丧失;而且修正是多余的——信赖域本身就引入了修正(解信赖域问题相当于分解 \(\nabla^2f+\lambda I\),\(\lambda\) 由半径决定)。结论:dogleg 最适合凸目标(Hessian 恒半正定);一般情形宜用下面的方法。

信赖域问题的精确解(PDF p.175–176)

用 Algorithm 4.4:重复分解 \(B_k+\lambda I\),经验上每次迭代平均解 1–3 个系统,代价不算过高。得到的算法很稳健——可期望收敛到极小点而不只是驻点。大规模时每次迭代解多个线性系统可能负担过重。

信赖域 Newton–CG(Trust-Region Newton–CG Method,PDF p.176–177)

在 Algorithm 4.3(CG–Steihaug)中取 \(B_k=\nabla^2f(x_k)\):对

\[B_kp_k=-\nabla f_k \tag{6.27}\]
做 CG,在 (i) 近似解超出半径,(ii) 达到所需精度,(iii) 遇负曲率(沿负曲率方向走到边界)时停止。是 Algorithm 6.1 的信赖域对应物。

  • 控制内层 CG 精度是降低代价的关键。在良态解附近信赖域约束不起作用,退化为 6.1 节的非精确牛顿,强制序列决定后期收敛速度。
  • 同样可 Hessian-free(自动微分或 (6.7))。
  • 优点:全局收敛(第一步沿 \(-\nabla f_k\) 即 Cauchy 点,之后 CG 只会改进模型值);无需矩阵分解,可利用稀疏性而不担心填充;CG 的核心是矩阵–向量乘,可并行;Hessian 正定时越接近解越逼近纯牛顿步,快速收敛。
  • 相比线搜索 Newton–CG:步长受信赖域控制;探索负曲率方向——经验上有益,有时能让迭代离开非极小驻点。
  • 局限:接受任何负曲率方向,即使模型下降微不足道。例(PDF 文本为 \(m(p)=10^{-3}p_1+10^{-4}p_1^2-p_2^2\),最速下降方向印作 \((10^{-3},0)^T\);按作者论述,自洽的写法应为 \(m(p)=10^{-3}p_1-10^{-4}p_1^2-p_2^2\),\(\|p\|\le1\),\(p=0\) 处最速下降方向 \((-10^{-3},0)^T\) 恰是负曲率方向):Algorithm 4.3 沿该方向走到边界,模型只降约 \(10^{-3}\);而沿 \(e_2\)(同为负曲率方向)可降 1。
  • 补救:Hessian 有负特征值时,方向应在最负特征值的特征向量上有显著分量,以快速离开非极小驻点。用 Lanczos 方法代替 CG:遇第一个负曲率方向后不终止,继续寻找“足够负”的曲率方向(更稳健但子问题更贵)。

Newton–CG 的预条件(Preconditioning the Newton–CG Method,PDF p.177–179)

  • Hessian 病态时无预条件 CG 可能低效甚至达不到精度,需预条件:找非奇异 \(D\) 使 \(D^{-T}B_kD^{-1}\) 的特征值分布更好。
  • 难点:预条件 CG 的迭代不再保持 \(\ell_2\) 范数单调增(定理 4.2 只对无预条件成立),不能一碰边界就停(后续迭代可能回到区域内)。
  • 解决:存在依赖于预条件子的加权范数使迭代单调增。考虑
    \[\min_p m_k(p)\quad\text{s.t. }\|Dp\|\le\Delta_k, \tag{6.28}\]
    令 \(\hat p=Dp\),\(\hat g_k=D^{-T}\nabla f_k\),\(\hat B_k=D^{-T}B_kD^{-1}\),化为标准形式 (6.24),对其直接用 CG;Algorithm 4.3 监控 \(\|\hat p\|=\|Dp\|\)(单调增),超过 \(\Delta_k\) 即终止。即用椭球信赖域与预条件子匹配。
  • 常用预条件子:不完全 Cholesky,\(B=LL^T-R\),\(L\) 的填充受限(如与 \(B\) 下三角同稀疏结构),\(R\) 为不精确部分。还需处理 Hessian 可能不定。

Algorithm 6.6(非精确修正 Cholesky):

(尺度化)T = diag(‖B e_i‖);B̄ ← T^{−1/2} B T^{−1/2};β ← ‖B̄‖
(计算平移以保证正定)if min_i b_ii > 0: α_0 ← 0  else α_0 ← β/2
for k = 0,1,2,...
  尝试不完全 Cholesky:L Lᵀ = B̄ + α_k I
  if 成功: 返回 L   else α_{k+1} ← max(2α_k, β/2)
end

取预条件子 \(D=L\)。MINPACK-2 的 NMTR 实现了此信赖域 Newton–CG;LANCELOT 也有使用略不同预条件的 Newton–CG。

信赖域牛顿法的局部收敛(Local Convergence,PDF p.179–182)

关键:证明接近解时信赖域约束最终不起作用,近似解在区域内部并越来越接近纯牛顿步。满足后者的步称渐近精确(asymptotically exact)。

定理 6.4:\(f\) 二阶 Lipschitz 连续可微,\(x_k\to x^*\)(满足二阶充分条件),充分大 \(k\) 时算法取 \(B_k=\nabla^2f(x_k)\),步 \(p_k\) 至少达到 Cauchy 下降(\(m_k(p_k)\le m_k(p_k^C)\)),且当 \(\|p_k^N\|\le\frac12\Delta_k\) 时渐近精确:

\[\|p_k-p_k^N\|=o(\|p_k^N\|). \tag{6.29}\]
则充分大 \(k\) 时信赖域约束不起作用。 证明:

  1. 无论 \(\|p_k^N\|\le\frac12\Delta_k\) 与否,都有 \(\|p_k\|\le2\|p_k^N\|\le2\|\nabla^2f_k^{-1}\|\|\nabla f_k\|\),即 \(\|\nabla f_k\|\ge\frac12\|p_k\|/\|\nabla^2f_k^{-1}\|\)。
  2. 由 Cauchy 下降估计 (4.34):\(m_k(0)-m_k(p_k)\ge c_1\|\nabla f_k\|\min(\Delta_k,\|\nabla f_k\|/\|B_k\|)\ge c_1\dfrac{\|p_k\|^2}{4\|\nabla^2f_k^{-1}\|^2\|\nabla^2f_k\|}\);由连续性,充分大 \(k\) 时 \(\ge c_3\|p_k\|^2\)(6.30),\(c_3=\dfrac{c_1}{8\|\nabla^2f(x^*)^{-1}\|^2\|\nabla^2f(x^*)\|}\)。
  3. 由 Hessian Lipschitz:\(|(f(x_k)-f(x_k+p_k))-(m_k(0)-m_k(p_k))|\le\frac L2\|p_k\|^3\),故
    \[|\rho_k-1|\le\frac{L}{2c_3}\|p_k\|\le\frac{L}{2c_3}\Delta_k. \tag{6.31}\]
  4. 半径只在 \(\rho_k<\frac14\) 时缩小,故存在阈值 \(\tilde\Delta\),\(\Delta_k\) 不会降到其下,\(\{\Delta_k\}\) 远离零;而 \(\|p_k^N\|\to0\),由 (6.29) \(\|p_k\|\to0\),约束最终不起作用。

引理 6.5:\(x_k\to x^*\)(满足二阶充分条件),用 \(B_k=\nabla^2f(x_k)\) 的 dogleg (6.25) 或二维子空间 (6.26):充分大 \(k\) 时模型无约束极小就是 \(p_k^N\),满足定理 6.4。Newton–CG 用 (6.3) 且 \(\eta_k\to0\)(加越界/负曲率终止)也满足。 证明:dogleg/子空间:\(p_k^N\) 在区域内、在 dogleg 路径上、在子空间内,且是区域内模型极小,故 \(p_k=p_k^N\)。Newton–CG:充分大 \(k\) 时 Hessian 正定,不会因负曲率停止;CG 迭代范数递增但不超过 \(\|p_k^N\|\),留在区域内;只能因 (6.3) 停止,于是 \(\|p_k-p_k^N\|\le\|\nabla^2f_k^{-1}\|\|r_k\|\le\eta_k\|\nabla^2f_k^{-1}\|\|\nabla f_k\|\le\eta_k\|\nabla^2f_k^{-1}\|\|\nabla^2f_k\|\|p_k^N\|\),满足 (6.29)。 Algorithm 4.4 也满足(牛顿步在区域内时 \(\lambda=0\)、\(p_k=p_k^N\) 最终满足任何合理终止准则)。

结论:最终取精确牛顿步的方法二次收敛;Newton–CG 的渐近速度由 \(\eta_k\to0\) 的速度决定(定理 6.2)。

注释与参考(PDF p.182)

迭代法求步长的牛顿法见 Sherman、Ortega–Rheinboldt;非精确牛顿讨论基于 Dembo–Eisenstat–Steihaug。修正 Cholesky 见 Gill–Murray–Wright、Dennis–Schnabel;Gershgorin 版见 Schnabel–Eskow;修正不定分解见 Cheng–Higham。另一种线搜索策略是计算负曲率方向并用于构造搜索方向(Moré–Sorensen、Goldfarb)。

习题概览(PDF p.182–183)

  • 6.1 编程:无线搜索的纯牛顿迭代,用 CG 求方向,选停止准则分别得线性/超线性/二次收敛;测试凸四次函数 \(f(x)=\frac12x^Tx+0.25\sigma(x^TAx)^2\)(6.32),\(A=\begin{bmatrix}5&1&0&0.5\\1&4&0.5&0\\0&0.5&3&0\\0.5&0&0&2\end{bmatrix}\),\(x_1=(\cos70°,\sin70°,\cos70°,\sin70°)^T\),\(\sigma=1\) 或更大(\(\sigma\) 控制偏离二次的程度)。
  • 6.2 证明线搜索 Newton–CG 方向总是下降方向。
  • 6.3 计算 (6.21) 两个 \(2\times2\) 块的特征值(各一正一负),验证 \(A\) 与 \(B\) 惯性相同。
  • 6.4 修正 Cholesky 对 \(\mathrm{diag}(-2,12,4)\) 的作用。
  • 6.5 编程:无预条件 CG–Steihaug 信赖域 Newton–CG,超线性参数,测试 (6.32);再构造在 \(x_0=0\) 附近有负曲率的函数观察负曲率步。
  • 6.6 证明定理 6.2;6.7 块对角矩阵的特征分解可逐块计算。

本章要点

  1. 非精确牛顿:\(\|r_k\|\le\eta_k\|\nabla f_k\|\);\(\eta_k\le\eta<1\) 线性,\(\eta_k\to0\) 超线性,\(\eta_k=O(\|\nabla f_k\|)\) 二次。
  2. 线搜索 Newton–CG(截断牛顿):CG 起点为 0,遇负曲率终止;Hessian-free(Hessian–向量积用差分或自动微分)。
  3. 修正牛顿法:\(B_k=\nabla^2f_k+E_k\);有界修正分解性质 ⇒ 全局收敛;修正最终消失 ⇒ 二次收敛。
  4. Hessian 修正:特征值修正(Frobenius 最优 / \(\tau I\) 平移)、加 \(\tau I\) 的 Cholesky 试探、修正 Cholesky(GMW,\(d_j=\max(|c_{jj}|,(\theta_j/\beta)^2,\delta)\))、Gershgorin、修正对称不定分解(Bunch–Parlett/Kaufman)。
  5. 信赖域牛顿:dogleg 适合凸问题;精确解法稳健;信赖域 Newton–CG 适合大规模且能利用负曲率,可用 Lanczos 增强;预条件需配合椭球信赖域。
  6. 局部收敛:渐近精确步使信赖域约束最终失效,恢复牛顿法的快速收敛。

与量化交易的关联

  • 协方差矩阵修复:样本相关/协方差矩阵因缺失数据成对估计、人工调整相关性或压力情景设定而不再半正定,是风险建模中的常见问题。(6.12) 的“特征值截断到 \(\delta\)”就是最常用的修复法(再重新标准化对角为 1,近似 Higham 最近相关矩阵思路);(6.13) 的整体平移 \(+\tau I\) 对应收缩到单位阵;修正 Cholesky 适合在需要 Cholesky 因子(蒙特卡洛生成相关正态随机数)时直接得到可用因子。
  • 极大似然估计与 Hessian:GARCH、Copula、随机波动率模型的 MLE 在远离最优时 Hessian 常不定;线搜索牛顿需修正 Hessian,否则可能走上升方向。估计标准误时用的是最优点处的 Hessian 逆,若最优点 Hessian 接近奇异(参数弱识别),本章“修正可能不消失、只线性收敛”的结论提示应检查模型可识别性。
  • 大规模 Hessian-free 方法:大规模组合优化(数千资产 + 非线性成本)、大规模 logistic 回归/因子选择等,可用截断牛顿(SciPy Newton-CG、trust-ncg),只需梯度与 Hessian–向量积,内存需求低。
  • 非精确求解的经济学:强制序列思想——远离最优时不必精确求解子问题——对日内高频再平衡的实时优化(时间预算有限)很有启发:先粗解保证下降,接近最优再提高精度。
  • 对称不定分解与 KKT 系统:带等式约束的组合优化 KKT 矩阵 \(\begin{bmatrix}\Sigma&A^T\\A&0\end{bmatrix}\) 是对称不定的,求解器(如内点法)内部用的正是 Bunch–Kaufman 型 \(LBL^T\) 分解,惯性检查用于判断是否为极小点(原书第 16 章展开)。

推荐习题

  • 6.1(强制序列与收敛速度的数值验证)。
  • 6.2(Newton–CG 方向的下降性证明,巩固第 5 章 CG 性质)。
  • 6.4(手算修正 Cholesky,理解对角修正行为)。
  • 6.5(信赖域 Newton–CG 与负曲率步的实验)。
  • 6.6(证明定理 6.2)。

第 7 章 导数计算(Calculating Derivatives,PDF p.184–211)

7.0 引言(PDF p.185–186)

多数非线性优化与非线性方程算法需要导数。手算可行时可让用户提供;函数太复杂时需自动计算或近似。三类方法:

  • 有限差分(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)。
  • 符号微分(symbolic differentiation):用 Mathematica、Maple、Macsyma 等对代数表达式符号运算。 本章讨论前两种。导数也用于最优后敏感性分析(post-optimal sensitivity analysis,最优解对参数/约束值小扰动的敏感度)以及非线性微分方程等。

7.1 有限差分导数近似(Finite-Difference Derivative Approximations,PDF p.186–196)

许多软件在用户不提供导数时自动做有限差分;结果虽近似,在很多情形已足够。

近似梯度(Approximating the Gradient,PDF p.186–189)

前向差分(forward / one-sided difference):

\[\frac{\partial f}{\partial x_i}(x)\approx\frac{f(x+\epsilon e_i)-f(x)}{\epsilon}, \tag{7.1}\]
整个梯度需 \(n+1\) 次函数值。依据:\(f(x+p)=f(x)+\nabla f(x)^Tp+\frac12p^T\nabla^2f(x+tp)p\)(7.2),若 \(\|\nabla^2f\|\le L\),则 \(|f(x+p)-f(x)-\nabla f(x)^Tp|\le(L/2)\|p\|^2\)(7.3);取 \(p=\epsilon e_i\):
\[\frac{\partial f}{\partial x_i}(x)=\frac{f(x+\epsilon e_i)-f(x)}{\epsilon}+\delta_\epsilon,\quad|\delta_\epsilon|\le(L/2)\epsilon. \tag{7.4}\]

步长选择与舍入误差:截断误差要求 \(\epsilon\) 小,但浮点运算有舍入误差。单位舍入 \(u\)(双精度约 \(10^{-16}\))是每次浮点运算相对误差的上界。粗略假设计算的 \(f\) 相对误差不超过 \(u\):\(|\mathrm{comp}(f(x))-f(x)|\le uL_f\)(\(L_f\) 为 \(|f|\) 的界),则总误差界

\[(L/2)\epsilon+2uL_f/\epsilon, \tag{7.5}\]
极小化得 \(\epsilon^2=4L_fu/L\)。问题尺度良好时 \(L_f/L\) 适中,取
\[\epsilon=\sqrt u \tag{7.6}\]
近乎最优(很多软件采用),总误差约 \(\sqrt u\)(约 \(10^{-8}\),即只有约一半有效数字)。

中心差分(central difference):

\[\frac{\partial f}{\partial x_i}(x)\approx\frac{f(x+\epsilon e_i)-f(x-\epsilon e_i)}{2\epsilon}, \tag{7.7}\]
需 \(2n\) 次函数值,约贵一倍。Hessian Lipschitz 时 \(f(x+p)=f(x)+\nabla f^Tp+\frac12p^T\nabla^2f(x)p+O(\|p\|^3)\)(7.8),取 \(p=\pm\epsilon e_i\) 相减得误差 \(O(\epsilon^2)\)。考虑求值误差后,最佳 \(\epsilon\approx u^{1/3}\)、精度约 \(u^{2/3}\)(习题 7.1;PDF 正文印作“\(\epsilon=u^{2/3}\)”,与习题 7.1 不一致,按习题为准)。多出的几位精度有时值得额外代价。

近似稀疏 Jacobian(Approximating a Sparse Jacobian,PDF p.189–193)

向量函数 \(r:\mathbb{R}^n\to\mathbb{R}^m\)(第 10 章残差、第 11 章非线性方程),Jacobian \(J(x)\) 见 (7.34)。由 Taylor:\(\|r(x+p)-r(x)-J(x)p\|\le(L/2)\|p\|^2\)(7.9)。

  • Jacobian–向量积(如非线性方程的非精确牛顿法所需):\(J(x)p\approx\dfrac{r(x+\epsilon p)-r(x)}{\epsilon}\)(7.10),精度 \(O(\epsilon)\);也可双侧。
  • 全 Jacobian 逐列:\(\dfrac{\partial r}{\partial x_i}(x)\approx\dfrac{r(x+\epsilon e_i)-r(x)}{\epsilon}\)(7.11),需 \(n+1\) 次 \(r\) 求值。稀疏时可同时估计多列,有时只需三四次。

例:

\[r(x)=\begin{bmatrix}2(x_2^3-x_1^2)\\3(x_2^3-x_1^2)+2(x_3^3-x_2^2)\\3(x_3^3-x_2^2)+2(x_4^3-x_3^2)\\\vdots\\3(x_n^3-x_{n-1}^2)\end{bmatrix} \tag{7.12}\]
每个分量只依赖两三个变量,\(n=6\) 时 Jacobian 三对角(7.13)。扰动 \(\epsilon e_1\) 只影响 \(r_1,r_2\);再加 \(\epsilon e_4\) 只影响 \(r_3,r_4,r_5\),两者互不干扰。取 \(p=\epsilon(e_1+e_4)\):\([r(x+p)]_{1,2}=[r(x+\epsilon e_1)]_{1,2}\)(7.14),\([r(x+p)]_{3,4,5}=[r(x+\epsilon e_4)]_{3,4,5}\)(7.15),于是
\[\begin{bmatrix}\partial r_1/\partial x_1\\\partial r_2/\partial x_1\end{bmatrix}\approx\frac{[r(x+p)-r(x)]_{1,2}}{\epsilon},\qquad\begin{bmatrix}\partial r_3/\partial x_4\\\partial r_4/\partial x_4\\\partial r_5/\partial x_4\end{bmatrix}\approx\frac{[r(x+p)-r(x)]_{3,4,5}}{\epsilon}. \tag{7.16–7.17}\]
(PDF 中 (7.17) 下标写作 \(\partial r_4/\partial x_3\) 等、并误称“Hessian 第四列”,按上下文应为 Jacobian 第 4 列的 \((3,4),(4,4),(5,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 incidence / intersection graph)\(\mathcal{G}\):\(n\) 个节点,若某个 \(r_j\) 同时依赖 \(x_i,x_k\)(即第 \(i,k\) 列在某行都有非零)则连边(图 7.1)。给节点着色,相邻节点不同色;每种颜色对应一个扰动向量 \(p=\epsilon(e_{i_1}+\cdots+e_{i_\ell})\)。最少着色难求,但有便宜的近似最优算法(Curtis–Powell–Reid,Coleman–Moré)。Newsam–Ramsdell:用更一般的扰动向量,至多 \(n_z\) 次额外求值(\(n_z\) 为每行最大非零数)。带状等结构已知最优着色;本例三色最优。

近似 Hessian(Approximating the Hessian,PDF p.193–194)

  • 有梯度无 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\) 个点;稀疏时跳过已知为零的元素。

近似稀疏 Hessian(Approximating a Sparse Hessian,PDF p.194–196)

Hessian 是 \(\nabla f\) 的 Jacobian,可直接用稀疏 Jacobian 技术,但忽略了对称性:估计了 \((i,j)\) 元也就有了 \((j,i)\) 元,利用对称性可大幅节省。 例:

\[f(x)=x_1\sum_{i=1}^ni^2x_i^2, \tag{7.22}\]
Hessian 呈“箭头形”(第一行、第一列与对角非零,7.23)。列关联图是完全图(第 1 行每列都非零),按 Jacobian 规则需 \(n+1\) 次梯度。利用对称:先用 \(p=\epsilon e_1\) 估第一列(即第一行);剩下只有对角元 \(2,\dots,6\),其图完全不连通,可同色,取
\[p=\epsilon(e_2+e_3+\cdots+e_6)=\epsilon(0,1,1,1,1,1)^T, \tag{7.24}\]
\(\nabla f\) 的第 \(i\) 分量只受 \(x_i\)(及 \(x_1\))扰动影响,故 \(\dfrac{\partial^2f}{\partial x_i^2}\approx\dfrac{\nabla f(x+\epsilon p)_i-\nabla f(x)_i}{\epsilon}\),\(i=2,\dots,6\)。总共只需在 \(x\) 及另外两点求梯度。 一般方法:用邻接图(adjacency graph:\(i\neq k\) 且 \(\partial^2f/\partial x_i\partial x_k\neq0\) 时连边),着色要求:相邻节点不同色,且任何长度为 3 的路径(\(i_1-i_2-i_3-i_4\))至少用三种颜色(Coleman–Moré)。

7.2 自动微分(Automatic Differentiation,PDF p.196–209)

利用函数的计算表示得到导数的解析值(非近似)。有的技术直接变换函数代码得到一般点的导数代码;有的记录特定点 \(x\) 处的计算过程再处理得导数。基础:任何函数都是一列一元/二元基本运算(加、乘、除、幂 \(a^b\);三角、指数、对数);链式法则

\[\nabla_xh(y(x))=\sum_{i=1}^m\frac{\partial h}{\partial y_i}\nabla y_i(x). \tag{7.25}\]
两种基本模式:前向与反向。

例(An Example,PDF p.197–198)

\[f(x)=(x_1x_2\sin x_3+e^{x_1x_2})/x_3. \tag{7.26}\]

计算图(图 7.2)引入中间变量:

\[x_4=x_1x_2,\ x_5=\sin x_3,\ x_6=e^{x_4},\ x_7=x_4x_5,\ x_8=x_6+x_7,\ x_9=x_8/x_3. \tag{7.27}\]
有向边 \(i\to j\) 时称 \(i\) 为 \(j\) 的父节点、\(j\) 为 \(i\) 的子节点;父节点值已知即可计算,计算从左到右流动,称前向扫描(forward sweep)。AD 工具自动识别中间量、构造计算图,用户无需手工分解。

前向模式(The Forward Mode,PDF p.198–199)

对每个中间变量同步计算沿方向 \(p\) 的方向导数:

\[D_px_i\overset{\text{def}}{=}(\nabla x_i)^Tp=\sum_{j=1}^3\frac{\partial x_i}{\partial x_j}p_j,\quad i=1,\dots,9. \tag{7.28}\]
目标 \(D_px_9=\nabla f(x)^Tp\)。自变量处 \(D_px_i=p_i\),\(p\) 称种子向量(seed vector)。例:\(x_7=x_4x_5\),
\[D_px_7=\frac{\partial x_7}{\partial x_4}D_px_4+\frac{\partial x_7}{\partial x_5}D_px_5=x_5D_px_4+x_4D_px_5. \tag{7.29}\]

  • 实现:不需同时存储所有节点的 \(x_i\)、\(D_px_i\)——节点的所有子节点算完后即可覆盖。软件给每个标量 \(w\) 关联 \(D_pw\),每次运算都做相应链式运算,如 \(z=w/y\):
    \[D_pz\leftarrow\frac1yD_pw-\frac{w}{y^2}D_py. \tag{7.30}\]
  • 全梯度:对 \(n\) 个种子 \(e_1,\dots,e_n\) 同时进行。代价可能显著:一次除法引起约 \(2n\) 次乘法和 \(n\) 次加法;存储可能增加 \(n\) 倍。计算早期很多 \(D_{e_j}x_i\) 为零,可用稀疏数据结构节省。
  • 实现方式:预编译器(把函数代码变换为扩展代码)或 C++ 等语言的运算符重载。

反向模式(The Reverse Mode,PDF p.199–202)

先完成 \(f\) 的计算,再反向扫描计算图,求 \(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{+}=\dfrac{\partial f}{\partial x_j}\dfrac{\partial x_j}{\partial x_i}\)(7.32)。所有子节点贡献完毕,\(\bar x_i=\partial f/\partial x_i\),节点“完成”(finalized),再向其父节点贡献。计算方向是从子到父,与求值相反。
  • 反向扫描用数值而非公式:前向扫描时不仅算 \(x_i\),还计算并存储每条边的偏导数值 \(\partial x_j/\partial x_i\)。

数值例:\(x=(1,2,\pi/2)^T\)(图 7.3,\(p(j,i)=\partial x_j/\partial x_i\)):\(x_4=2\),\(x_5=1\),\(x_6=e^2\),\(x_7=2\),\(x_8=2+e^2\),\(x_9=(4+2e^2)/\pi\);边值 \(p(4,1)=2\),\(p(4,2)=1\),\(p(5,3)=\cos(\pi/2)=0\),\(p(6,4)=e^2\),\(p(7,4)=x_5=1\),\(p(7,5)=x_4=2\),\(p(8,6)=p(8,7)=1\),\(p(9,8)=1/x_3=2/\pi\),\(p(9,3)=-x_8/x_3^2=-(2+e^2)/(\pi/2)^2\)。 反向:\(\bar x_9=1\);

\[\bar x_3\mathrel{+}=\frac{\partial f}{\partial x_9}\frac{\partial x_9}{\partial x_3}=-\frac{2+e^2}{(\pi/2)^2}=\frac{-8-4e^2}{\pi^2},\qquad\bar x_8\mathrel{+}=\frac{1}{\pi/2}=\frac2\pi. \tag{7.33}\]
节点 3 还等待子节点 5 的贡献;节点 8 只有子节点 9,完成;然后 \(\bar x_6\mathrel{+}=2/\pi\),\(\bar x_7\mathrel{+}=2/\pi\),进而更新节点 4、5,最终
\[\begin{bmatrix}\bar x_1\\\bar x_2\\\bar x_3\end{bmatrix}=\nabla f(x)=\begin{bmatrix}(4+4e^2)/\pi\\(2+2e^2)/\pi\\(-8-4e^2)/\pi^2\end{bmatrix}.\]
(验证:\(\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\)。)

  • 主要优点:对标量函数 \(f:\mathbb{R}^n\to\mathbb{R}\),梯度的额外算术量至多为函数求值的 4–5 倍,与 \(n\) 无关(例:(7.33a) 需两次乘、一次除、一次加,(7.33b) 一次除一次加,约为前向中那一次除法的 5 倍)。前向模式可能需 \(n\) 倍。对向量函数 \(r:\mathbb{R}^n\to\mathbb{R}^m\),随 \(m\) 增大两者代价趋近。
  • 缺点:需存储整个计算图(每次基本运算存一个节点:中间结果、指向一两个父节点的指针、边上偏导数);反向时按写入逆序读取,访问模式简单;可通过运算符重载实现(ADOL-C),反向扫描作为一次函数调用。但存储可能巨大:每节点 20 字节,在 100 megaflop 机器上求值 1 秒的函数,图可达 2 GB。检查点(checkpointing):对图的片段做部分前向/反向扫描,按需重算而非全部存储,以额外算术换存储(Griewank 等;Odyssee)。

向量函数与部分可分性(Vector Functions and Partial Separability,PDF p.203–204)

  • 非线性最小二乘与非线性方程中 \(r:\mathbb{R}^n\to\mathbb{R}^m\),计算图最右列有 \(m\) 个无子节点,Jacobian
    \[J(x)=\left[\frac{\partial r_j}{\partial x_i}\right]_{j=1..m,\ i=1..n}. \tag{7.34}\]
  • 部分可分(partially separable)函数:
    \[f(x)=\sum_{i=1}^{n_e}f_i(x), \tag{7.35}\]
    每个元素函数(element function)\(f_i\) 只依赖 \(x\) 的少数分量。令 \(r(x)=(f_1(x),\dots,f_{n_e}(x))^T\),则
    \[\nabla f(x)=J(x)^Te,\quad e=(1,\dots,1)^T. \tag{7.36}\]
    \(J\) 多数列非零元少,可用图着色高效计算再恢复梯度。Gay 开发了自动检测部分可分分解的方法,用户无需额外信息即可利用其效率(第 9 章有对应拟牛顿方法)。
  • 约束优化中,同时计算 \(f\) 与约束 \(c_i\)(\(r(x)=(f(x),[c_j(x)]_{j\in\mathcal{I}\cup\mathcal{E}})\))可共享公共子表达式(如图 7.2 中 \(x_4\) 被 \(x_6,x_7\) 共享),减少总工作量。

计算向量函数的 Jacobian(PDF p.204–205)

  • 前向模式:给种子 \(p\),最右节点得 \(D_pr_j=(\nabla r_j)^Tp\),组装即 \(J(x)p\)(Jacobian–向量积)。全 Jacobian 取 \(p=e_1,\dots,e_n\);稀疏时用着色选种子。算术增加倍数约等于种子数。
  • 反向模式:选种子 \(q\in\mathbb{R}^m\),对标量函数 \(r(x)^Tq\) 做反向:\(\nabla[r(x)^Tq]=J(x)^Tq\)(Jacobian 转置–向量积)。在 \(m\) 个因变量节点的 \(\bar x\) 上置 \(q_1,\dots,q_m\),结束时自变量节点含 \(J(x)^Tq\) 的分量。全 Jacobian 取 \(q=e_1,\dots,e_m\);稀疏时对 \(J^T\) 着色。算术增加不超过种子数的 5 倍;存储不超过标量情形。
  • 两模式可组合:前向种子揭示部分列,反向种子揭示其余行。
  • 非精确牛顿等只需反复算 \(J(x)p\),一次前向扫描即可,代价与函数求值相当。

计算 Hessian:前向模式(PDF p.205–207)

对种子对 \((p,q)\) 定义

\[D_{pq}x_i=p^T(\nabla^2x_i)q, \tag{7.37}\]
前向扫描中与 \(x_i\)、\(D_px_i\) 同步计算;自变量处 \(D_{pq}=0\);最右节点得 \(p^T\nabla^2f(x)q\)。传播规则:\(x_i=x_j+x_k\) 时
\[D_px_i=D_px_j+D_px_k,\quad D_{pq}x_i=D_{pq}x_j+D_{pq}x_k; \tag{7.38}\]
一元运算 \(x_i=L(x_j)\):
\[x_i=L(x_j),\quad D_px_i=L'(x_j)D_px_j,\quad D_{pq}x_i=L''(x_j)(D_px_j)(D_qx_j)+L'(x_j)D_{pq}x_j. \tag{7.39}\]
需同时累积 \(D_p\)、\(D_q\)。

  • 稠密 Hessian:所有单位向量对 \((e_j,e_k)\),\(k\le j\),共 \(n(n+1)/2\) 对(对称只算下三角);已知稀疏结构时只算可能非零的位置。算术增加约为 \(1+n+N_z(\nabla^2f)\) 的小倍数,每个累积量每节点一个存储位(可覆盖)。
  • 只需 Hessian–向量积 \(\nabla^2f(x)q\):算 \(D_{e_1}x_i,\dots,D_{e_n}x_i\)、\(D_qx_i\) 及 \(D_{e_1q}x_i,\dots,D_{e_nq}x_i\),最后节点给出 \(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\) 的一、二阶导;可推广到更高阶导(Bischof–Corliss–Griewank)。

计算 Hessian:反向模式(PDF p.207–208)

求 \(\nabla^2f(x)q\):前向模式同时算 \(f\) 与 \(\nabla f(x)^Tq\)(累积 \(x_i\)、\(D_qx_i\)),再对计算出的函数 \(\nabla f(x)^Tq\) 做标准反向扫描,自变量节点得 \(\dfrac{\partial}{\partial x_i}(\nabla f(x)^Tq)=[\nabla^2f(x)q]_i\)。算术增加与 \(n\) 无关:前向部分为小倍数,反向再乘至多 5,总计约为 \(f\) 的 12 倍。全 Hessian 取 \(q=e_1,\dots,e_n\),至多 \(12n\) 倍;稀疏已知结构时用着色,约 \(12N_c(\nabla^2f)\) 倍(\(N_c\) 为种子数)。

当前局限(Current Limitations,PDF p.208–209)

  • 若 \(f\) 的求值依赖 PDE 数值解,计算值 \(\hat f(x)=f(x)+\tau(x)\) 含截断误差;\(|\tau|\) 小但 \(\tau'\) 未必小,导数误差可能很大(有限差分同样受影响);用分段有理函数近似三角函数时也有类似问题。
  • 代码中的分支:病态例 \(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 是日益成熟的增强技术,使优化算法能用于复杂函数,也便于解释最优解(敏感性);但不是万灵药,用户仍需思考导数计算。

注释与参考(PDF p.209)

AD 工具:ODYSSEE、ADIFOR、ADIC、ADOL-C;发展趋势是前向/反向的“混合模式”(Bischof–Haghighat)。综述:Griewank–Corliss 主编文集(至 1991)、Berz–Bischof–Corliss–Griewank 续集(含 Corliss–Rall 及详尽文献目录)、Griewank 关于 AD 在优化中应用的综述。部分可分函数梯度计算见 Bischof 等;Hessian 计算见 Gay 等。稀疏 Hessian 高效估计:Coleman–Moré(图着色语言)之前已有 Powell–Toint 的高效方案;稀疏 Hessian/Jacobian 估计软件见 Coleman–Garbow–Moré。

习题概览(PDF p.209–211)

  • 7.1 中心差分最优 \(\epsilon=u^{1/3}\),可达精度约 \(u^{2/3}\)。
  • 7.2 推导 Hessian–向量积的中心差分;7.3 验证 (7.21)。
  • 7.4 Hessian 对角非零时,邻接图是 \(\nabla f\) 关联图的子图;7.5 画 (7.22) 的邻接图,验证“节点 1 一色、其余一色”有效;7.6 给定稀疏结构找三色方案。
  • 7.7 跟踪 (7.26) 的前向模式(题中指标与变量数与正文不完全一致);7.8 推导加法、指数、正切、幂 \(s^t\) 的前向规则;7.9 验证图 7.3 的边值并完成反向扫描、给出节点完成顺序;7.10 乘法与余弦的反向规则及工作量比较;7.11 减、乘、除的 \(D_p\)、\(D_{pq}\) 规则;7.12 验证 (7.39)。
  • 7.13 \(f(x)=\frac12[x^Tx+(a^Tx)^2]\),统计计算 \(f\)、\(\nabla f\)、\(\nabla^2f\)、\(\nabla^2f(x)p\) 的运算量。

本章要点

  1. 前向差分误差 \(O(\epsilon)\),最佳 \(\epsilon\approx\sqrt u\),精度约 \(\sqrt u\);中心差分误差 \(O(\epsilon^2)\),最佳 \(\epsilon\approx u^{1/3}\),精度约 \(u^{2/3}\),代价加倍。
  2. Jacobian–向量积、Hessian–向量积只需一次额外函数/梯度求值。
  3. 稀疏 Jacobian/Hessian:图着色组合扰动向量,大幅减少求值次数;Hessian 还可利用对称性(邻接图 + 长度 3 路径三色规则)。
  4. 自动微分:前向模式(方向导数随求值传播,全梯度代价约 \(n\) 倍);反向模式(伴随变量反向累积,梯度代价至多 4–5 倍且与 \(n\) 无关,但需存储计算图,可用检查点)。
  5. 向量函数:前向得 \(Jp\),反向得 \(J^Tq\);部分可分函数 \(\nabla f=J^Te\)。
  6. Hessian:前向二阶传播 \(D_{pq}\);前向+反向得 \(\nabla^2f\,q\) 约 12 倍代价。
  7. 局限:截断误差的导数、代码分支。

与量化交易的关联

  • 定价敏感度(Greeks):有限差分是计算 Delta、Gamma、Vega 的最常用方法;本章的步长分析直接适用——前向差分步长取 \(\sqrt u\cdot\max(|x|,1)\) 量级;蒙特卡洛定价时函数值本身带统计噪声(远大于 \(u\)),最优步长要按噪声水平而非机器精度重新推导((7.5) 中把 \(u L_f\) 换成噪声标准差),且必须使用公共随机数,否则差分毫无意义。
  • AAD(伴随算法微分):反向模式在衍生品定价中称 AAD,可在约 4–5 倍单次定价成本下得到对上百个市场参数(收益率曲线节点、波动率曲面格点)的全部敏感度,与参数个数无关——这是 XVA 和大型衍生品账簿风险计算的核心技术,本章给出其原理。检查点技术对应长路径蒙特卡洛的内存问题。
  • 分支陷阱:Payoff 中的指示函数(数字期权、障碍期权)不连续,AD 和有限差分都会出错,需光滑化(smoothing)或似然比方法——对应本章“当前局限”中的分支例子。
  • 组合优化与风险归因:现代框架(PyTorch/JAX)的反向模式自动微分可直接对复杂组合目标(含非线性成本、风险预算)求梯度;风险贡献 \(w_i(\Sigma w)_i/\sigma_p\) 本质就是组合波动率对权重的梯度(欧拉分解),AD 可自动得到。
  • 稀疏 Jacobian 着色:利率曲线自举/多曲线校准中,每个校准工具只依赖少数曲线节点,Jacobian 带状稀疏,着色技术可显著减少重定价次数。

推荐习题

  • 7.1(中心差分最优步长推导,直接用于 Greeks 计算的步长选择)。
  • 7.2、7.3(Hessian–向量积与纯函数值 Hessian 差分公式)。
  • 7.5、7.6(稀疏 Hessian 的图着色)。
  • 7.9、7.10(手工执行反向模式,理解 AAD)。
  • 7.13(运算量统计,比较各种导数计算方式的代价)。

第 8 章 拟牛顿法(Quasi-Newton Methods,PDF p.212–219,本块只覆盖至 8.1 节“BFGS 方法的性质”开头)(续见下一块)

8.0 引言与历史(PDF p.213–214)

  • 1950 年代中期,Argonne 国家实验室的物理学家 W.C. Davidon 用坐标下降法做长时间优化计算,当时计算机不稳定,总在算完前崩溃。他为加速迭代发明了第一个拟牛顿算法——非线性优化中最具革命性的思想之一。Fletcher 与 Powell 很快证明它比现有方法快且可靠得多,一夜之间改变了非线性优化。此后二十年出现大量变体和数百篇论文。讽刺的是,Davidon 的论文当年未被接受发表,作为技术报告存在三十多年,直到 1991 年刊于 SIAM Journal on Optimization 创刊号。
  • 拟牛顿法与最速下降一样每步只需梯度;通过测量梯度变化构造足以产生超线性收敛的模型。相对最速下降改进巨大,尤其在难题上;因不需二阶导,有时比牛顿法更高效。现在软件库包含多种拟牛顿算法(无约束、约束、大规模)。本章讨论中小规模,第 9 章讨论大规模推广。
  • 自动微分削弱了拟牛顿法的吸引力(消除手算二阶导的繁琐与出错风险),但只是有限程度;拟牛顿法在许多问题上仍有竞争力。

8.1 BFGS 方法(The BFGS Method,PDF p.214–219,未完)

最流行的拟牛顿算法,以发现者 Broyden、Fletcher、Goldfarb、Shanno 命名。本节推导 BFGS 及其近亲 DFP。

模型与迭代:当前点二次模型

\[m_k(p)=f_k+\nabla f_k^Tp+\tfrac12p^TB_kp, \tag{8.1}\]
\(B_k\) 为 \(n\times n\) 对称正定、每步更新;模型在 \(p=0\) 的值和梯度与 \(f_k\)、\(\nabla f_k\) 一致。凸二次模型的极小点
\[p_k=-B_k^{-1}\nabla f_k \tag{8.2}\]
作为搜索方向,
\[x_{k+1}=x_k+\alpha_kp_k, \tag{8.3}\]
\(\alpha_k\) 满足 Wolfe 条件 (3.6)。与线搜索牛顿法类似,区别在于用近似 Hessian \(B_k\)。

割线方程的推导:Davidon 提出不每次重算 \(B_k\),而是简单地更新它以反映最近一步测得的曲率。新模型 \(m_{k+1}(p)=f_{k+1}+\nabla f_{k+1}^Tp+\frac12p^TB_{k+1}p\);合理要求:\(m_{k+1}\) 的梯度在最近两个迭代点 \(x_k,x_{k+1}\) 处与 \(f\) 的梯度一致。在 \(x_{k+1}\) 处自动满足;在 \(x_k\) 处:\(\nabla m_{k+1}(-\alpha_kp_k)=\nabla f_{k+1}-\alpha_kB_{k+1}p_k=\nabla f_k\),整理得

\[B_{k+1}\alpha_kp_k=\nabla f_{k+1}-\nabla f_k. \tag{8.4}\]
记
\[s_k=x_{k+1}-x_k,\quad y_k=\nabla f_{k+1}-\nabla f_k, \tag{8.5}\]
得割线方程(secant equation)
\[B_{k+1}s_k=y_k. \tag{8.6}\]

曲率条件:对称正定 \(B_{k+1}\) 把 \(s_k\) 映为 \(y_k\),只有在

\[s_k^Ty_k>0 \tag{8.7}\]
时可能((8.6) 左乘 \(s_k^T\) 即见)。\(f\) 强凸时对任意两点成立(习题);非凸时不总成立,需通过线搜索强制。Wolfe 或强 Wolfe 条件保证 (8.7):由 (3.6b),\(\nabla f_{k+1}^Ts_k\ge c_2\nabla f_k^Ts_k\),故
\[y_k^Ts_k\ge(c_2-1)\alpha_k\nabla f_k^Tp_k, \tag{8.8}\]
因 \(c_2<1\) 且 \(p_k\) 下降,右端为正。(这是拟牛顿法必须用 Wolfe 而非仅 Armijo 线搜索的根本原因。)

唯一性——最近矩阵问题:曲率条件满足时割线方程有无穷多解:对称矩阵有 \(n(n+1)/2\) 个自由度,割线方程只给 \(n\) 个条件;正定性给 \(n\) 个不等式(所有主子式为正),仍不能吸收剩余自由度。要求 \(B_{k+1}\) 在某种意义上最接近 \(B_k\):

\[\min_B\|B-B_k\|\quad\text{s.t. }B=B^T,\ Bs_k=y_k. \tag{8.9}\]
不同范数给出不同拟牛顿法。采用便于求解且导致尺度不变方法的加权 Frobenius 范数
\[\|A\|_W\equiv\|W^{1/2}AW^{1/2}\|_F,\quad \|C\|_F^2=\sum_{i,j}c_{ij}^2, \tag{8.10}\]
权 \(W\) 为任一满足 \(Wy_k=s_k\) 的矩阵,具体可取 \(W=\bar G_k^{-1}\),\(\bar G_k\) 为平均 Hessian
\[\bar G_k=\int_0^1\nabla^2f(x_k+\tau\alpha_kp_k)d\tau, \tag{8.11}\]
由 Taylor 定理 \(y_k=\bar G_k\alpha_kp_k=\bar G_ks_k\)(8.12)。此权使范数无量纲(解不依赖问题单位)。

DFP 公式:(8.9) 的唯一解为

\[B_{k+1}=(I-\gamma_ky_ks_k^T)B_k(I-\gamma_ks_ky_k^T)+\gamma_ky_ky_k^T,\quad\gamma_k=\frac1{y_k^Ts_k}. \tag{8.13}\]
由 Davidon 1959 年提出,Fletcher 与 Powell 研究、实现并推广,故称 DFP。逆近似 \(H_k=B_k^{-1}\) 便于用矩阵–向量乘计算方向;由 Sherman–Morrison–Woodbury 公式得
\[H_{k+1}=H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{y_k^Ts_k}. \tag{8.14}\]
右端后两项为秩一矩阵,\(H_k\) 做秩二修正;(8.13) 也是 \(B_k\) 的秩二修正。拟牛顿更新的基本思想:不每步从头重算,而用简单修正把最新观测到的目标信息与当前近似中已有的知识结合。

BFGS 公式:DFP 很有效但很快被 BFGS 取代,后者被认为是所有拟牛顿公式中最有效的。推导:对逆 \(H_k\) 施加类似条件——\(H_{k+1}\) 对称正定且满足 \(H_{k+1}y_k=s_k\),

\[\min_H\|H-H_k\|\quad\text{s.t. }H=H^T,\ Hy_k=s_k, \tag{8.15}\]
加权 Frobenius 范数的权 \(W\) 满足 \(Ws_k=y_k\)(如 \(W=\bar G_k\))。唯一解:
\[H_{k+1}=(I-\rho_ks_ky_k^T)H_k(I-\rho_ky_ks_k^T)+\rho_ks_ks_k^T, \tag{8.16}\]
\[\rho_k=\frac1{y_k^Ts_k}. \tag{8.17}\]
(DFP 与 BFGS 互为“对偶”:交换 \(B\leftrightarrow H\)、\(s\leftrightarrow y\)。)

初始矩阵 \(H_0\):没有万能公式。可利用问题信息(如 \(x_0\) 处有限差分近似 Hessian 的逆),或取单位阵、单位阵的倍数(倍数反映变量尺度)。

Algorithm 8.1(BFGS 方法):

给定 x_0,容差 ε > 0,逆 Hessian 近似 H_0;k ← 0
while ‖∇f_k‖ > ε
  p_k = −H_k ∇f_k                                     (8.18)
  x_{k+1} = x_k + α_k p_k,α_k 由满足 Wolfe 条件 (3.6) 的线搜索给出
  s_k = x_{k+1} − x_k,y_k = ∇f_{k+1} − ∇f_k
  按 (8.16) 计算 H_{k+1}
  k ← k+1
end
  • 每步 \(O(n^2)\) 运算(加函数与梯度求值),无解线性方程组或矩阵–矩阵乘等 \(O(n^3)\) 运算。稳健,超线性收敛,对多数实用目的已足够快。牛顿法虽二次收敛,但每步要解线性方程组,代价更高;BFGS 更重要的优势是不需二阶导。
  • \(B\) 形式:对 (8.16) 用 Sherman–Morrison–Woodbury 得
    \[B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}. \tag{8.19}\]
    朴素实现需解 \(B_kp_k=-\nabla f_k\),代价 \(O(n^3)\),不适合无约束极小化;但通过更新 \(B_k\) 的 Cholesky 因子可得更便宜的实现(后文讨论)。

BFGS 方法的性质(Properties of the BFGS Method,PDF p.219,本块读到此处)

数值例:Rosenbrock 函数 (2.23),初值 \((-1.2,1)\),三种方法都用 Wolfe 步长,将梯度范数降到 \(10^{-5}\):最速下降 5264 次迭代,BFGS 34 次,非精确牛顿 21 次。最后几步的 \(\|x_k-x^*\|\):

最速下降 BFGS 牛顿
1.827e-04 1.70e-03 3.48e-02
1.826e-04 1.17e-03 1.44e-02
1.824e-04 1.34e-04 1.82e-04
1.823e-04 1.01e-06 1.17e-08

可见最速下降线性且极慢(比率接近 1),BFGS 超线性,牛顿二次。

正定性的保持:(8.15) 并未显式要求正定,但 \(H_k\) 正定 ⇒ \(H_{k+1}\) 正定。证明:由 (8.8),\(y_k^Ts_k>0\),(8.16)(8.17) 有定义。对任意 \(z\neq0\),

\[z^TH_{k+1}z=w^TH_kw+\rho_k(z^Ts_k)^2\ge0,\quad w=z-\rho_ky_k(s_k^Tz).\]
右端为零仅当 \(s_k^Tz=0\),但此时 \(w=z\neq0\),第一项为正。故 \(H_{k+1}\) 正定。

尺度不变性:要使更新公式对变量变换不变,(8.9a)(8.15a) 的目标也须不变,权矩阵 \(W\) 的选取保证了这一点。其他 \(W\) 给出其他公式,但尽管搜寻很多,没有找到显著优于 BFGS 的公式。

BFGS 用于二次函数时有很多有趣性质,将在更一般的 Broyden 族(BFGS 是其特例)背景下讨论。(续见下一块:8.1 节实现、8.2 SR1、8.3 Broyden 族、8.4 收敛分析等。)

本块第 8 章已读部分要点

  1. 拟牛顿模型 (8.1)–(8.3),割线方程 \(B_{k+1}s_k=y_k\),曲率条件 \(s_k^Ty_k>0\) 由 Wolfe 线搜索保证。
  2. 最近矩阵问题 + 加权 Frobenius 范数(权为平均 Hessian)⇒ DFP(对 \(B\))与 BFGS(对 \(H\))唯一解,均为秩二更新。
  3. BFGS 逆形式每步 \(O(n^2)\),无需二阶导,超线性收敛;\(H_k\) 正定性自动保持。
  4. Rosenbrock 上:最速下降 5264 步、BFGS 34 步、牛顿 21 步。

与量化交易的关联(第 8 章已读部分)

  • 最常用的通用优化器:SciPy minimize 默认(无约束时)就是 BFGS;L-BFGS-B(第 9 章)是带权重界约束的组合优化、MLE 估计的首选。理解“必须用 Wolfe 线搜索以保证 \(s_k^Ty_k>0\)”能解释为何自写 BFGS 配简单回溯时会出现 Hessian 近似失去正定、方向变成上升方向的问题。
  • MLE 估计与标准误:GARCH 等模型常用 BFGS 估计,并直接拿 BFGS 末步的 \(H_k\) 作为信息矩阵逆来算标准误——这是常见误用:\(H_k\) 只在搜索方向上逼近真 Hessian 逆(Dennis–Moré 条件),并不保证整体收敛到 \(\nabla^2f(x^*)^{-1}\),标准误应用数值或解析 Hessian 重新计算。
  • 初始矩阵与尺度:参数量级差异大(如 GARCH 中 \(\omega\sim10^{-6}\)、\(\alpha,\beta\sim0.1\))时,\(H_0\) 取单位阵效果差,应按变量尺度设对角 \(H_0\) 或先做参数变换。

推荐习题(第 8 章)

本章习题位于下一块(原书 p.220–221,PDF 约 p.240–241),由下一块负责概括。可先自行完成:验证强凸函数对任意两点满足 \(s^Ty>0\);由 (8.16) 用 Sherman–Morrison–Woodbury 推导 (8.19)。


全块说明

本块(PDF p.1–219)覆盖:前置部分;第 1–7 章全文(含各章习题);第 8 章开头至 8.1 节“BFGS 方法的性质”首段(原书 p.199)。第 8 章其余内容(8.1 实现、8.2 SR1、8.3 Broyden 族、8.4 收敛分析、注释与习题)在下一块。