第 03 章 线搜索方法
学习目标
读完本章,你应当能够:
- 写出 Armijo 充分下降条件、曲率条件、Wolfe 与强 Wolfe 条件、Goldstein 条件,解释每个条件排除了什么样的"坏步长",并会用回溯法实现充分下降。
- 陈述并证明 Zoutendijk 定理,用它说明"方向不与梯度垂直 + Wolfe 步长 ⇒ 全局收敛"。
- 定量说明最速下降的线性收敛率 \(\left(\frac{\kappa-1}{\kappa+1}\right)^2\) 及其在病态问题上的后果。
- 理解 Dennis–Moré 条件为何是拟牛顿法超线性收敛的充要条件,掌握牛顿法局部二次收敛的证明。
- 知道坐标下降法的优缺点,以及实用线搜索(插值、初始步长、Algorithm 3.2 + zoom)的设计要点。
- 在 SciPy 中排查 "line search failed",并解释参数估计中最速下降为何常常"跑不动"。
读前导读
这一章在解决什么问题
Solver 每一轮迭代都要做两个决定:往哪个方向改权重,改多少。第 2 章讲了"往哪个方向"(最速下降、牛顿、拟牛顿);本章讲"改多少",也就是步长。这件事听起来琐碎,实际上决定了算法是稳稳收敛、原地打转,还是一步冲过头。SciPy 报 "line search failed"、Solver 报"无法找到可行解"或在某个点附近来回跳,很多时候都是步长这一环出了问题。
可以用调仓来想象。你已经确定"应该减配股票、加配债券"这个方向,接下来要决定调 1% 还是 10%。调太少,每次只改善一点点,可能永远到不了目标;调太多,可能越过最优点,反而更差。本章的 Wolfe 条件就是两条"验收标准":第一条(Armijo)要求"实际改善至少达到预测改善的一小部分",防止冲过头;第二条(曲率条件)要求"走到的位置坡度已经明显变缓",防止走得太短。
本章第二个主题是收敛性理论:在什么条件下,算法保证最终会停在梯度为零的点上(Zoutendijk 定理),以及停下来要多久(最速下降为什么慢、牛顿为什么快)。量化实战里你会看到一个很现实的后果:用最速下降估计 \(t\) 分布的自由度,跑完 2 万步、不报错,但给出的厚尾参数是错的。
需要先想起来的数学
- 一元函数 \(\phi(\alpha)\) 与链式法则:把多元函数限制在一条直线上,\(\phi(\alpha)=f(x_k+\alpha p_k)\) 就是普通的一元函数。它的导数 \(\phi'(\alpha)=\nabla f(x_k+\alpha p_k)^Tp_k\),用的是链式法则:外层对 \(x\) 求梯度,内层 \(x\) 对 \(\alpha\) 求导得 \(p_k\)。\(\phi'(0)<0\) 就是"出发时沿这个方向在下坡"。见 第 00 册第 02 章 导数与泰勒展开。
- 级数收敛:\(\sum_k a_k<\infty\)(非负项)意味着 \(a_k\to 0\)。这就像一笔永续年金现值有限,后期现金流必须越来越小。Zoutendijk 定理正是用这个逻辑推出梯度趋于零。见 第 00 册第 04 章 级数与收敛。
- Lipschitz 连续:\(\|\nabla f(x)-\nabla f(\tilde x)\|\le L\|x-\tilde x\|\),意思是"梯度变化不会比距离变化快 \(L\) 倍",即函数的弯曲程度有上限。一元时相当于 \(|f''|\le L\)。例如 \(f=x^2\) 的导数 \(2x\) 满足 \(L=2\)。见 第 00 册第 01 章 函数极限与连续。
- 特征值与条件数:二次函数 \(\frac12x^TQx\) 的等高线是椭圆,椭圆各轴的"陡峭程度"就是 \(Q\) 的特征值 \(\lambda_1\le\dots\le\lambda_n\)。条件数 \(\kappa=\lambda_n/\lambda_1\) 越大,椭圆越扁。见 第 00 册第 06 章 线性代数速成。
- \(\liminf\) 与 \(o(\cdot)\):\(\liminf_{k\to\infty}a_k=0\) 读作"下极限为零",意思是"有无穷多项能任意接近 0",但不要求全部项都趋于 0。\(o(\|p\|)\) 是"比 \(\|p\|\) 高阶的小量"。见 第 00 册第 07 章 概率中的分析工具。
怎么读这一章
必读:3.1、3.2.1(Wolfe 条件,配合原书图 3.3–3.5 看)、3.2.3(回溯法,最常用)、3.3 的 Zoutendijk 定理陈述与推论、3.4.1 的结论(不必推 Kantorovich 不等式)。3.5 的 Dennis–Moré 条件读懂"含义"即可;牛顿法二次收敛的证明值得第二遍细读,后面几章反复用到同样的技巧。3.6 坐标下降读"量化视角"一段。3.7 的插值公式和 Algorithm 3.2/3.3 属于实现细节,第一次可以跳过,但请读最后的"排查 line search failed"。
3.1 线搜索框架
线搜索方法的每一步是
其中 \(p_k\) 是搜索方向,\(\alpha_k>0\) 是步长(step length)。成功与否取决于两者都选得好。多数方法要求 \(p_k\) 是下降方向(\(p_k^T\nabla f_k<0\)),并且取
\(B_k\) 是对称非奇异矩阵:最速下降取 \(B_k=I\),牛顿法取 \(B_k=\nabla^2f(x_k)\),拟牛顿法中 \(B_k\) 由低秩公式逐步更新。若 \(B_k\) 正定,则 \(p_k^T\nabla f_k=-\nabla f_k^TB_k^{-1}\nabla f_k<0\),必为下降方向。
本章回答两个问题:步长怎么选(3.2–3.3 节、3.7 节),以及这样选了之后算法收敛吗、收敛多快(3.4–3.6 节)。
3.2 步长应满足的条件
选步长有一个基本权衡:希望 \(\alpha_k\) 使 \(f\) 大幅下降,又不愿意花太多时间去选。理想的选择是一元函数
的全局极小点,但这太贵了——即使只求局部极小到中等精度,也要很多次函数和梯度求值。实用策略是非精确线搜索(inexact line search):以最小代价得到"足够的"下降。典型过程分两阶段:包围阶段(bracketing)找到一个包含合适步长的区间,二分或插值阶段在区间内算出一个好步长。
只要求函数值下降是不够的。原书图 3.2 给了一个例子:函数的极小值是 \(-1\),但若迭代得到的函数值序列是 \(5/k\),它每步都在下降,却收敛到 0 而不是 \(-1\)。每一步下降的量太小,加起来也到不了最优值。所以需要"充分下降"。
3.2.1 Wolfe 条件
充分下降条件(sufficient decrease condition,又称 Armijo 条件):
它要求下降量与步长 \(\alpha\) 和方向导数 \(\nabla f_k^Tp_k\) 成比例。右端记为 \(l(\alpha)\),是一条斜率为 \(c_1\nabla f_k^Tp_k<0\) 的直线。因为 \(c_1<1\),在 \(\alpha\) 很小时这条线在 \(\phi\) 的图像之上,条件即 \(\phi(\alpha)\le l(\alpha)\)(原书图 3.3)。实践中 \(c_1\) 取得很小,如 \(10^{-4}\)。
白话解释:\(\alpha\nabla f_k^Tp_k\) 是一阶泰勒展开"预测"的下降量(斜率 × 步长)。Armijo 条件说:实际下降至少要达到预测下降的 \(c_1\) 倍。数值例:\(\phi'(0)=\nabla f_k^Tp_k=-2\),试 \(\alpha=1\),预测下降 2,\(c_1=10^{-4}\) 时只要求实际下降 \(\ge 0.0002\)。这个门槛很低,几乎只排除"走过头、函数值反而上升或几乎没降"的步长。 这和预算执行的检查类似:你不要求实际利润达到预测值,只要求它不能离谱地低于预测——否则说明预测模型(线性近似)在这么远的地方已经不可信了,需要把步子收小。 但这个条件单独使用不够:所有足够小的 \(\alpha\) 都满足它,算法可能因步长太短而停滞。为排除过短步长,加上曲率条件(curvature condition):
左端正是 \(\phi'(\alpha_k)\),条件要求 \(\phi'(\alpha_k)\ge c_2\phi'(0)\)。直观理解:如果在 \(\alpha_k\) 处斜率仍然很负,说明沿这个方向继续走还能显著下降,不应停;如果斜率只是略负甚至为正,说明这个方向上已没多少下降空间,可以停了(图 3.4)。典型取值:牛顿/拟牛顿法 \(c_2=0.9\);非线性共轭梯度法 \(c_2=0.1\)(要求更接近一维极小)。
两者合称 Wolfe 条件:
注意满足 Wolfe 条件的步长不一定接近 \(\phi\) 的极小点(图 3.5)。若希望排除离驻点很远的点,用强 Wolfe 条件(strong Wolfe conditions):
区别在于不再允许 \(\phi'(\alpha_k)\) 是很大的正数。
白话解释:把两条 Wolfe 条件放在一起,用一元函数 \(\phi(\alpha)\) 的斜率来理解。出发点斜率 \(\phi'(0)\) 是负的(下坡)。取 \(c_2=0.9\)、\(\phi'(0)=-2\): 曲率条件 (3.6b) 要求 \(\phi'(\alpha)\ge -1.8\),即"到达处的下坡坡度至少比出发时缓 10%"。如果还是 \(-1.95\) 这么陡,说明走得太短,继续走还能大幅下降。 强 Wolfe (3.7b) 要求 \(|\phi'(\alpha)|\le 1.8\),即 \(-1.8\le\phi'(\alpha)\le 1.8\),额外排除了"已经越过谷底、在对面山坡上很陡地往上爬"的点。\(c_2\) 越小,这个区间越窄,步长越接近一维极小点。 为什么用"斜率变缓"而不是"函数值"来判断步长太短?因为在很短的步长上,函数值几乎没变,看不出问题;而斜率仍然很陡,一眼就能看出"还没走完"。 引理 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(\alpha)\) 有下界,而直线 \(l(\alpha)=f(x_k)+\alpha c_1\nabla f_k^Tp_k\) 斜率为负、无下界,所以二者至少相交一次。设最小的交点为 \(\alpha'>0\):
所有小于 \(\alpha'\) 的步长都满足充分下降。由中值定理,存在 \(\alpha''\in(0,\alpha')\) 使
联立 (3.8)(3.9):
(因为 \(c_1<c_2\) 且 \(\nabla f_k^Tp_k<0\))。所以 \(\alpha''\) 严格满足 Wolfe 两个条件,由光滑性其附近有一个区间也满足。又因 (3.10) 左端为负,\(|\phi'(\alpha'')|=c_1|\phi'(0)|<c_2|\phi'(0)|\),强 Wolfe 条件也在该区间成立。\(\square\)
Wolfe 条件在广义上是尺度不变的:目标函数乘一个正常数,或对变量做仿射变换,都不改变哪些步长满足条件。它适用于多数线搜索方法,对拟牛顿法尤其重要——第 8 章会看到,曲率条件保证 \(s_k^Ty_k>0\),这正是 BFGS 更新保持正定所需要的。
3.2.2 Goldstein 条件
右侧不等式就是充分下降;左侧从下方限制步长,排除过短的步(图 3.6)。缺点是左侧不等式可能把 \(\phi\) 的所有极小点都排除掉。Goldstein 条件与 Wolfe 条件的收敛理论相似,常用于牛顿型方法,但不适合拟牛顿法——它不能保证 \(s_k^Ty_k>0\),无法维持 Hessian 近似的正定性。
3.2.3 充分下降与回溯
只用充分下降条件也行,前提是候选步长按回溯(backtracking)方式产生:
Procedure 3.1(回溯线搜索)
- 选 \(\bar\alpha>0\),\(\rho\in(0,1)\),\(c\in(0,1)\),令 \(\alpha\leftarrow\bar\alpha\);
- 当 \(f(x_k+\alpha p_k)>f(x_k)+c\alpha\nabla f_k^Tp_k\) 时,令 \(\alpha\leftarrow\rho\alpha\);
- 返回 \(\alpha_k=\alpha\)。
要点:
- 牛顿法和拟牛顿法取 \(\bar\alpha=1\);最速下降和共轭梯度法可取其他初值。
- 有限次试探后必然终止:\(\alpha\) 足够小时充分下降一定成立(因为 \(p_k\) 是下降方向)。
- 收缩因子 \(\rho\) 可以每次不同(比如由带保护的插值决定),只要保持在 \([\rho_{lo},\rho_{hi}]\subset(0,1)\) 内。
- 回溯为什么能避免步长过短?被接受的 \(\alpha_k\) 要么就是初值 \(\bar\alpha\),要么与上一个被拒绝(太长)的 \(\alpha_k/\rho\) 只差一个常数因子——所以它不会比"必要的"短太多。
- 回溯很适合牛顿法(第 6 章),不太适合拟牛顿法和共轭梯度法(它们需要曲率条件)。
3.3 线搜索方法的收敛性
全局收敛不仅需要好的步长,还需要好的方向。关键量是 \(p_k\) 与最速下降方向 \(-\nabla f_k\) 的夹角 \(\theta_k\):
定理 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 连续:
则
证明:由曲率条件 (3.6b) 两边减去 \(\nabla f_k^Tp_k\):
由 Lipschitz 条件:\((\nabla f_{k+1}-\nabla f_k)^Tp_k\le\alpha_kL\|p_k\|^2\)。两式合并得步长下界
代入充分下降条件 (3.6a):
从 0 到 \(k\) 累加:
\(f\) 有下界,所以级数收敛。\(\square\)
推导拆解:把证明拆成四步,并补上省略的中间式。 第 1 步(曲率条件变形):(3.6b) 是 \(\nabla f_{k+1}^Tp_k\ge c_2\nabla f_k^Tp_k\),两边同减 \(\nabla f_k^Tp_k\) 得 \((\nabla f_{k+1}-\nabla f_k)^Tp_k\ge(c_2-1)\nabla f_k^Tp_k\)。右边是"负数乘负数",为正:走完这一步后,沿 \(p_k\) 的斜率至少回升了这么多。 第 2 步(Lipschitz 给上界):\(x_{k+1}-x_k=\alpha_kp_k\),用柯西–施瓦茨不等式 \(u^Tv\le\|u\|\|v\|\) 和 (3.13):\((\nabla f_{k+1}-\nabla f_k)^Tp_k\le\|\nabla f_{k+1}-\nabla f_k\|\|p_k\|\le L\alpha_k\|p_k\|^2\)。斜率回升量受"弯曲程度上限 × 走的距离"约束。 第 3 步(合并得步长下界):\(L\alpha_k\|p_k\|^2\ge(c_2-1)\nabla f_k^Tp_k\),两边除以 \(L\|p_k\|^2\)。含义:要让斜率回升到位,必须走够一定距离,所以步长不会太短。 第 4 步(代入 Armijo):(3.6a) 给出 \(f_{k+1}\le f_k+c_1\alpha_k\nabla f_k^Tp_k\)。因为 \(\nabla f_k^Tp_k<0\),把 \(\alpha_k\) 换成它的下界会让右边更大,不等式仍成立:\(f_{k+1}\le f_k-c_1\frac{1-c_2}{L}\frac{(\nabla f_k^Tp_k)^2}{\|p_k\|^2}\)。再用 (3.12) 的定义,\((\nabla f_k^Tp_k)^2/\|p_k\|^2=\cos^2\theta_k\|\nabla f_k\|^2\)。 最后,每一步的下降量相加,类似把每期现金流加总:\(f\) 从 \(f_0\) 出发、永远不低于某个下界,所以所有下降量之和有限。 这个证明值得记住它的结构:曲率条件 + Lipschitz ⇒ 步长不会太短;充分下降 ⇒ 每步下降量至少正比于 \(\cos^2\theta_k\|\nabla f_k\|^2\);函数有下界 ⇒ 总下降量有限。 后面分析信赖域、共轭梯度时会反复看到同样的思路。假设条件也不苛刻:\(f\) 无下界时问题本身无意义;梯度 Lipschitz 一般被局部收敛理论中的光滑性假设所蕴含。用 Goldstein 或强 Wolfe 条件也有同样的结论。
推论 Zoutendijk 条件 (3.14) 蕴含
如果方向与梯度的夹角始终远离 \(90°\),即存在 \(\delta>0\) 使
那么
特别地,用 Wolfe 或 Goldstein 线搜索的最速下降法,梯度收敛到零(此时 \(\cos\theta_k=1\))。
白话解释:(3.14) 说一个非负数列的和有限,所以每一项 \(\cos^2\theta_k\|\nabla f_k\|^2\) 必须趋于 0。这个乘积趋于 0 有两种可能:梯度变小(我们想要的),或者 \(\cos\theta_k\to0\),即搜索方向越来越接近与梯度垂直(沿等高线走,几乎不下坡)。条件 (3.17) 排除了第二种可能,于是只剩第一种。 \(\theta_k\) 是搜索方向与"最陡下坡方向"的夹角。\(\cos\theta_k=1\) 是正对下坡,\(\cos\theta_k=0\) 是横着走。只要方向始终"有一定的下坡分量",Wolfe 步长就能保证最终走到平地。 本书说一个算法全局收敛(globally convergent),就是指满足 (3.18)。这是一般线搜索方法能得到的最强结果:它只保证迭代被驻点吸引,不保证收敛到极小点。要加强为收敛到局部极小,必须利用 Hessian 的负曲率信息(第 4 章的信赖域方法能做到)。
牛顿型方法 若 \(B_k\) 正定且条件数一致有界,\(\|B_k\|\|B_k^{-1}\|\le M\),则(练习 3)
从而 \(\lim\|\nabla f_k\|=0\)(3.20)。也就是说,\(B_k\) 正定、条件数有界、步长满足 Wolfe 条件,牛顿法和拟牛顿法就是全局收敛的。第 6 章的"有界修正分解"正是为满足这个条件而设计的。
较弱的结果 对共轭梯度法等,只能证明
即只有一个子列的梯度趋于零。证明套路是反证:若 \(\|\nabla f_k\|\ge\gamma>0\)(3.22),由 (3.16) 必有 \(\cos\theta_k\to0\)(3.23);只要能证明存在一个子列使 \(\cos\theta_{k_j}\) 远离零,就得到矛盾。第 5 章用这个套路证明 Fletcher–Reeves 方法的收敛性。
一个实用推论:任何算法,只要 (i) 每步都使 \(f\) 下降,(ii) 每 \(m\) 步做一次满足 Wolfe/Goldstein 条件的最速下降步,就满足 (3.20)。其余 \(m-1\) 步可以做"更聪明"的事;偶尔插入的最速下降步也许进展不大,但它是全局收敛的保险。
3.4 收敛速度(一):最速下降
看起来,只要保证 \(p_k\) 不与梯度趋于垂直(或定期走最速下降步)就万事大吉。但如果每步检查 \(\cos\theta_k\),一旦小于 \(\delta\) 就把方向往最速下降拉(所谓角度检验),会带来两个问题:(1) Hessian 病态时,好的方向本来就可能几乎与梯度垂直,\(\delta\) 选得不当会阻碍快速收敛;(2) 破坏拟牛顿法的不变性。快速收敛与全局收敛常常相互冲突:最速下降是全局收敛的典范,但很慢;纯牛顿法在解附近极快,但远离解时步长甚至可能不是下降方向。算法设计的挑战在于兼得两者。
3.4.1 二次函数上的精确分析
考虑理想情形——二次目标 + 精确线搜索:
梯度 \(\nabla f(x)=Qx-b\),唯一极小点 \(x^*\) 是 \(Qx=b\) 的解。沿 \(-\nabla f_k\) 方向的精确步长为(练习 4)
迭代为
原书图 3.7 展示了它的典型行为:等高线是轴沿 \(Q\) 特征向量的椭圆,迭代呈之字形(zigzag)前进。
引入加权范数 \(\|x\|_Q^2=x^TQx\),利用 \(Qx^*=b\) 可以验证
所以 \(Q\)-范数误差就是函数值误差。进一步可推出(习题 3.7)
定理 3.3 精确线搜索的最速下降法用于强凸二次函数 (3.24) 时,
其中 \(0<\lambda_1\le\cdots\le\lambda_n\) 是 \(Q\) 的特征值。
证明只需在 (3.28) 中用 Kantorovich 不等式(习题 3.8):对任意 \(x\ne0\),
而 \(1-\frac{4\lambda_n\lambda_1}{(\lambda_n+\lambda_1)^2}=\left(\frac{\lambda_n-\lambda_1}{\lambda_n+\lambda_1}\right)^2\)。
用条件数 \(\kappa(Q)=\lambda_n/\lambda_1\) 表示,收敛因子是 \(\left(\frac{\kappa-1}{\kappa+1}\right)^2\)。所以:
- 函数值线性收敛。所有特征值相等(等高线是圆)时一步收敛。
- \(\kappa\) 越大,椭圆越扁,之字形越严重,收敛越慢。(3.29) 是最坏情形的界,但原书指出 \(n>2\) 时它能相当准确地反映实际行为(Akaike 的概率分析表明,对大多数初始点最坏界都近似成立)。
推导拆解:精确步长 (3.25) 怎么来?令 \(g=\nabla f_k\),\(\phi(\alpha)=f(x_k-\alpha g)\)。展开二次函数得 \(\phi(\alpha)=f_k-\alpha g^Tg+\tfrac12\alpha^2g^TQg\)(一次项用了 \(\nabla f_k=Qx_k-b\),二次项来自 \(\tfrac12(\alpha g)^TQ(\alpha g)\))。对 \(\alpha\) 求导令其为零:\(-g^Tg+\alpha g^TQg=0\),解出 \(\alpha=g^Tg/g^TQg\)。 \(\|x\|_Q=\sqrt{x^TQx}\) 是"按 \(Q\) 加权的长度",就像用协方差矩阵衡量组合偏离的跟踪误差 \(\sqrt{d^T\Sigma d}\),而不是简单的权重差之和。(3.27) 说这个加权长度的平方的一半,正好等于函数值离最优值还差多少。
金融直觉:把收敛因子 \(\left(\frac{\kappa-1}{\kappa+1}\right)^2\) 代几个数感受一下:\(\kappa=2\) 时约 0.11,每步误差缩到 11%;\(\kappa=100\) 时约 0.96;\(\kappa=10^4\) 时约 0.9996,要约 5700 步误差才缩小 10 倍。在组合优化里,\(\kappa\) 大往往是因为有两只高度相关的股票:它们的"多空价差"方向方差很小(特征值小),"同涨同跌"方向方差很大(特征值大)。最速下降在这种"长条形峡谷"里会反复横跳,迟迟找不准价差头寸。
3.4.2 一般非线性函数
定理 3.4 设 \(f\) 二阶连续可微,精确线搜索的最速下降迭代收敛到 \(x^*\),且 \(\nabla^2f(x^*)\) 正定。则
(渐近成立),其中 \(\lambda_i\) 是 \(\nabla^2f(x^*)\) 的特征值。非精确线搜索一般不会改善这个速度。
原书的数值例:\(\kappa=800\),\(f(x_1)=1\),\(f(x^*)=0\),原书称 1000 次迭代后函数值仍约为 0.08。
关于这个数字的说明:若严格按 (3.29) 的平方因子计算,\(\left(\frac{799}{801}\right)^{2\times1000}\approx e^{-5}\approx0.0067\);原书的 0.08 对应的是未平方的因子 \(\left(\frac{799}{801}\right)^{1000}\approx e^{-2.5}\approx0.082\)。两种算法得到的数值不同,但结论不变:即使 Hessian 的条件不算太坏,最速下降也可能慢到不可接受——上千步之后函数值误差仍在 \(10^{-2}\) 到 \(10^{-3}\) 量级。本章量化实战会对此做数值验证。
量化中协方差矩阵的条件数常在 \(10^3\)–\(10^5\) 量级(高度相关的资产、近似共线的因子、波动率差异巨大的资产类别),用最速下降求解组合或风险模型问题会非常慢,应改用牛顿/拟牛顿法或做预条件(第 5 章)。
3.5 收敛速度(二):拟牛顿法与牛顿法
3.5.1 拟牛顿法与 Dennis–Moré 条件
拟牛顿方向为
\(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^2f(x^*)\) 正定,且
则 (i) 存在 \(k_0\),当 \(k>k_0\) 时单位步长 \(\alpha_k=1\) 满足 Wolfe 条件;(ii) 若 \(k>k_0\) 时都取 \(\alpha_k=1\),则 \(\{x_k\}\) 超线性收敛。
两点说明:
- 为什么要 \(c_1\le\frac12\)?对二次函数,精确极小点处 \(\phi(\alpha^*)=\phi(0)+\frac12\alpha^*\phi'(0)\),若 \(c_1>\frac12\),充分下降条件会把这个极小点排除掉,单位步长可能永远不被接受。
- 对拟牛顿方向 (3.30),条件 (3.31) 等价于
这是一个令人惊喜的结论:并不需要 \(B_k\) 收敛到真实 Hessian,只要 \(B_k\) 沿搜索方向 \(p_k\) 越来越准,就能超线性收敛。
白话解释:(3.31) 的分子 \(\nabla f_k+\nabla^2f_kp_k\) 是什么?牛顿步正是让它等于零的那个 \(p\)(见第 2 章 (2.14) 的推导)。所以这个量衡量"你的步离牛顿步差多远",除以 \(\|p_k\|\) 变成相对误差。Dennis–Moré 条件就是:相对误差趋于零,即步子越来越像牛顿步。 用风险模型打个比方:你不需要整个协方差矩阵都估准,只需要在你实际要调仓的那个方向上把风险估准,调仓决策就和用真实协方差做出的决策几乎一样。拟牛顿法的割线方程恰好就是在"刚走过的方向"上校准 \(B_k\)。 定理 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) 成立。
证明思路:记牛顿步 \(p_k^N=-\nabla^2f_k^{-1}\nabla f_k\)。先证明 (3.32) 等价于
事实上
而 \(\|\nabla^2f_k^{-1}\|\) 在 \(x^*\) 附近有界,\(\nabla^2f_k\to\nabla^2f(x^*)\),所以 (3.32) ⇒ (3.33);反方向两边乘 \(\nabla^2f_k\) 即得。然后利用下文牛顿法的二次收敛估计 (3.37):
再由 \(\|p_k\|=O(\|x_k-x^*\|)\) 得 \(\|x_{k+1}-x^*\|=o(\|x_k-x^*\|)\),即超线性收敛。反方向类似。\(\square\)
第 8 章会证明 BFGS 等拟牛顿方法在合理条件下确实满足 (3.32)。
3.5.2 牛顿法
牛顿方向
Hessian 未必正定,\(p_k^N\) 未必是下降方向,所以本章许多结论不能直接套用。第 6 章给出两种全局化途径:线搜索(必要时修正 Hessian)与信赖域。这里只讨论局部收敛:在 \(\nabla^2f(x^*)\) 正定的解附近,Hessian 也正定;若步长最终总取 1,则二次收敛。
定理 3.7 设 \(f\) 二阶可微,Hessian 在满足二阶充分条件的解 \(x^*\) 的邻域内 Lipschitz 连续(常数 \(L\))。考虑 \(x_{k+1}=x_k+p_k^N\)。则:
- 若初始点足够接近 \(x^*\),迭代收敛到 \(x^*\);
- 收敛速度是二次的;
- 梯度范数 \(\|\nabla f_k\|\) 二次收敛到零。
证明:利用 \(\nabla f(x^*)=0\) 和牛顿步的定义,
由 (2.5),\(\nabla f_k-\nabla f_*=\int_0^1\nabla^2f(x_k+t(x^*-x_k))(x_k-x^*)dt\),于是
\(x_k\) 充分接近 \(x^*\) 时 \(\|\nabla^2f_k^{-1}\|\le2\|\nabla^2f(x^*)^{-1}\|\),于是
推导拆解:(3.35) 怎么得到?从 \(x_k+p_k^N-x^*\) 出发,代入 \(p_k^N=-\nabla^2f_k^{-1}\nabla f_k\),并在前两项前面"乘一个 \(\nabla^2f_k^{-1}\nabla^2f_k=I\)":\(x_k-x^*-\nabla^2f_k^{-1}\nabla f_k=\nabla^2f_k^{-1}[\nabla^2f_k(x_k-x^*)-\nabla f_k]\)。最后因为 \(\nabla f_*=\nabla f(x^*)=0\),把 \(\nabla f_k\) 写成 \(\nabla f_k-\nabla f_*\),这是为了下一步能用积分形式 (2.5)。 (3.36) 的含义:方括号里是"用 \(x_k\) 处的 Hessian 预测的梯度变化"减"实际梯度变化"。实际变化是沿线段上各点 Hessian 的平均,预测只用了起点的 Hessian。两者之差取决于 Hessian 沿线段变化了多少:Lipschitz 条件说变化量 \(\le L\cdot t\|x_k-x^*\|\),对 \(t\) 从 0 到 1 积分得到 \(\tfrac12\),再乘一个 \(\|x_k-x^*\|\),于是出现平方。 这就是二次收敛的来源:牛顿法的误差只来自"Hessian 在这一步内的变化",而这个变化本身与步长成正比,所以新误差 ∝ 旧误差的平方。类比:用久期 + 凸性估债券价格,误差来自"凸性本身随收益率变化",是 \((\Delta y)^3\) 量级。
只要初值满足 \(\tilde L\|x_0-x^*\|<1\),归纳即得收敛与二次速度。对梯度,利用 \(\nabla f_k+\nabla^2f_kp_k^N=0\):
牛顿方向使 (3.31) 中的比值恒为零,所以由定理 3.5,当 \(k\) 充分大时 Wolfe(以及 Goldstein)条件接受单位步长。结论:总是先试单位步长的牛顿法实现,最终都会取 \(\alpha_k=1\),从而获得局部二次收敛。 这也是计量软件做极大似然估计时,最后几步总能"瞬间"收敛的原因。
3.6 坐标下降法
坐标下降法(coordinate descent)依次以坐标方向 \(e_1,\dots,e_n\) 为搜索方向:固定其余变量,对一个变量极小化(或至少降低)\(f\),\(n\) 步后再从头循环。也叫交替变量法(method of alternating variables)。
- 简单直观,但实践中可能很低效:原书图 3.8 显示在一个二维二次函数上,几步之后横竖移动都进展甚微。
- 带精确线搜索的坐标下降可以无限迭代而永远不接近梯度为零的点(Powell 的例子),对比之下最速下降保证 \(\|\nabla f_k\|\to0\)。更一般地,沿任意一组固定的线性无关方向循环搜索都不保证全局收敛。原因是梯度可能越来越接近与坐标方向垂直,\(\cos\theta_k\) 足够快地趋于零,Zoutendijk 条件成立而 \(\nabla f_k\) 不趋于零。
- 即使收敛,通常也比最速下降慢得多,变量越多差距越大。
- 但它仍然有用:(1) 不需要计算梯度;(2) 变量弱耦合(loosely coupled)时速度可以接受。
- 变体:来回扫描 \(e_1,\dots,e_n,e_{n-1},\dots,e_1,\dots\);或在一轮坐标步后沿首末点连线再搜索一次(Hooke–Jeeves 等模式搜索方法基于此思想),其中部分变体是全局收敛的。
量化视角:Lasso / Elastic Net 因子选择的标准算法(glmnet)就是坐标下降。它之所以在这里很好用,一是 \(\ell_1\) 惩罚是可分的、每个坐标子问题有闭式解(软阈值);二是因子标准化后相关性不高时,变量正是"弱耦合"的。反过来,当候选因子高度相关时,glmnet 会明显变慢——这正是本节的结论。风险平价的循环坐标算法也属于这一类。
3.7 步长选择算法
本节讨论怎样实际算出满足上述条件的步长。记 \(\phi(\alpha)=f(x_k+\alpha p_k)\)(3.38),假设 \(\phi'(0)<0\),只在 \(\alpha>0\) 上搜索。
若 \(f\) 是凸二次函数 \(f(x)=\frac12x^TQx+b^Tx+c\),沿射线的一维极小点有解析式
一般非线性函数则需要迭代。线搜索对所有非线性优化方法的稳健性和效率都有很大影响。只用函数值的方法可能低效;利用梯度可以直接判断是否已找到合适步长(Wolfe/Goldstein 条件)。好消息是:接近解时,第一个试探步往往就满足条件。
所有过程都分两阶段:包围阶段找到一个区间 \([a,b]\);选择阶段(zoom)缩小区间并用插值猜测极小点位置。
3.7.1 插值
把 Procedure 3.1 增强一下:生成一列递减的 \(\alpha_i\),每个不比前一个小太多,直到满足充分下降
设计目标是尽量少算导数。
- 若初始猜测 \(\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}\]
- 若 \(\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}\)。
- 必要时重复,每次用 \(\phi(0)\)、\(\phi'(0)\) 和最近两个函数值做三次插值。
- 保护措施(safeguard):若新的 \(\alpha_i\) 离前一个太近,或比前一个小太多,就重置为 \(\alpha_i=\alpha_{i-1}/2\),以保证每步有合理进展、最终步长不会太小。
推导拆解:(3.41) 的思路是"用已知的三条信息拼一条抛物线"。设 \(\phi_q(\alpha)=a\alpha^2+b\alpha+c\)。\(\phi_q(0)=\phi(0)\) 给出 \(c=\phi(0)\);\(\phi_q'(0)=\phi'(0)\) 给出 \(b=\phi'(0)\);\(\phi_q(\alpha_0)=\phi(\alpha_0)\) 解出 \(a=[\phi(\alpha_0)-\phi(0)-\phi'(0)\alpha_0]/\alpha_0^2\)。抛物线极小点在 \(-b/(2a)\),即 (3.42)。 分母里的 \(\phi(\alpha_0)-\phi(0)-\phi'(0)\alpha_0\) 是"实际函数值"减"直线预测值"。\(\alpha_0\) 被 Armijo 拒绝说明实际值远高于直线预测,所以这个量为正,抛物线开口向上,极小点存在。 数值例:\(\phi(0)=0\),\(\phi'(0)=-1\),\(\alpha_0=1\),\(\phi(1)=0.5\)。则 \(a=1.5\),\(\alpha_1=1/3\)。它利用了"实际值比直线预测(\(-1\))高出 1.5、弯曲很明显"这条信息,而简单折半只会机械地试 0.5。
上述策略假设导数比函数值贵得多。实际上,方向导数常常可以几乎免费地与函数值一起算出(第 7 章的自动微分),这时可用最近两个点的 \(\phi\) 和 \(\phi'\) 做三次插值(对曲率变化明显的函数是很好的模型):
三次插值能使步长序列二次收敛到一维极小点。
3.7.2 初始步长
- 牛顿与拟牛顿法:总是先试 \(\alpha_0=1\),确保在条件允许时取单位步,发挥快速收敛性。
- 最速下降、共轭梯度等方向尺度不好的方法,需要利用当前信息猜测:
- 假设本步的一阶变化与上一步相同:\(\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)\) 做二次插值:
\[\alpha_0=\frac{2(f_k-f_{k-1})}{\phi'(0)}.\tag{3.44}\]若迭代超线性收敛,这个比值趋于 1。取 \(\alpha_0\leftarrow\min(1,1.01\alpha_0)\),单位步最终总会被试到并接受,从而在牛顿/拟牛顿法中也能观察到超线性收敛。
3.7.3 满足强 Wolfe 条件的线搜索算法
下面的算法对任意 \(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
逻辑依据:若 (i) \(\alpha_i\) 违反充分下降,或 (ii) \(\phi(\alpha_i)\ge\phi(\alpha_{i-1})\),或 (iii) \(\phi'(\alpha_i)\ge0\),那么区间 \((\alpha_{i-1},\alpha_i)\) 内一定含有满足强 Wolfe 条件的步长。试探步长单调增加,但传给 zoom 的两个参数顺序可能不同。外推选 \(\alpha_{i+1}\) 可以用插值或取常数倍,关键是增长要足够快,以便有限步内达到 \(\alpha_{\max}\)。
zoom\((\alpha_{lo},\alpha_{hi})\) 维护三个不变式:(a) \(\alpha_{lo}\) 与 \(\alpha_{hi}\) 之间含有满足强 Wolfe 条件的步长;(b) \(\alpha_{lo}\) 是迄今满足充分下降的步长中函数值最小者;(c) \(\phi'(\alpha_{lo})(\alpha_{hi}-\alpha_{lo})<0\)(从 \(\alpha_{lo}\) 朝 \(\alpha_{hi}\) 走函数值先下降)。
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}\)。
实现要点
- 插值步要有保护,避免过于靠近区间端点。
- 接近解时 \(f(x_k)\) 与 \(f(x_{k-1})\) 在有限精度下可能无法区分,要加入停止测试:例如若干次(典型为 10 次)试探都无法降低函数值就停止。
- 完整、稳健的实现很难写,原书建议使用公开的高质量代码,尤其是 Moré–Thuente 的实现。SciPy 的
scipy.optimize.line_search(以及minimize中 BFGS、CG、L-BFGS-B 的线搜索)就是这一类强 Wolfe 线搜索。 - 强 Wolfe 与普通 Wolfe 的代价:对"松"的线搜索(\(c_1=10^{-4}\),\(c_2=0.9\))两者工作量相近。强 Wolfe 的优势是减小 \(c_2\) 就能直接控制搜索质量(迫使 \(\alpha\) 更接近一维极小),这对最速下降和非线性共轭梯度很重要,所以实现强 Wolfe 的例程适用面更广。
排查 "line search failed":SciPy 报这个错误时,最常见的原因是 (1) 梯度写错了,\(p_k\) 实际不是下降方向;(2) 目标函数有噪声或不光滑(例如回测收益、带随机数的模拟),在小步长下函数值无法可靠下降;(3) 已经非常接近解,函数值的变化淹没在舍入误差中。第 (1) 条可用第 7 章的有限差分做梯度检查。
3.8 量化实战:线搜索在参数估计中的表现
我们做三个实验:
- A(原书习题 3.1):回溯线搜索的最速下降与牛顿法求解 Rosenbrock 函数,初值取 \((1.2,1.2)\) 和较难的 \((-1.2,1)\)。
- B:验证定理 3.3 的收敛率,并核对原书 \(\kappa=800\) 的数值例。
- C:用学生 \(t\) 分布的极大似然估计拟合厚尾日收益,参数 \(\theta=(\mu,\log\sigma,\log\nu)\)(取对数以去掉正值约束,变成无约束问题),比较"BFGS + Wolfe 线搜索"与"最速下降 + Wolfe 线搜索"。
import numpy as np
from scipy.optimize import line_search
from scipy.special import gammaln, digamma
# ---------- Part A: 回溯线搜索 + 最速下降 / 牛顿法,Rosenbrock(原书习题 3.1) ----------
def rosen(x): return 100 * (x[1] - x[0]**2)**2 + (1 - x[0])**2
def rosen_g(x): return np.array([-400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0]), 200 * (x[1] - x[0]**2)])
def rosen_h(x): return np.array([[1200 * x[0]**2 - 400 * x[1] + 2, -400 * x[0]], [-400 * x[0], 200.0]])
def backtracking(f, x, fx, g, p, alpha=1.0, rho=0.5, c=1e-4):
"""Procedure 3.1:缩短步长直到满足充分下降(Armijo)条件"""
while f(x + alpha * p) > fx + c * alpha * (g @ p):
alpha *= rho
return alpha
def run(direction, x0, tol=1e-8, maxit=50000):
x = np.array(x0, float); steps = []
for k in range(maxit):
g = rosen_g(x)
if np.linalg.norm(g) < tol: return k, x, steps
p = direction(x, g)
a = backtracking(rosen, x, rosen(x), g, p)
steps.append(a); x = x + a * p
return maxit, x, steps
sd = lambda x, g: -g
def newton(x, g):
H = rosen_h(x)
try:
np.linalg.cholesky(H); return -np.linalg.solve(H, g)
except np.linalg.LinAlgError: # Hessian 非正定时退回最速下降(第 6 章有更好的修正)
return -g
for x0 in [(1.2, 1.2), (-1.2, 1.0)]:
k_sd, x_sd, _ = run(sd, x0)
k_nt, x_nt, st = run(newton, x0)
print("x0=%-12s 最速下降 %5d 步 -> %s | 牛顿 %2d 步 -> %s,最后 5 步步长 %s" %
(x0, k_sd, np.round(x_sd, 6), k_nt, np.round(x_nt, 6), st[-5:]))
# ---------- Part B: 最速下降收敛率 (3.29) 的验证与原书 κ=800 的数值例 ----------
rng = np.random.default_rng(0)
n, kappa = 50, 800.0
lam = np.geomspace(1, kappa, n); Qm, _ = np.linalg.qr(rng.standard_normal((n, n)))
Q = Qm @ np.diag(lam) @ Qm.T; b = rng.standard_normal(n); xs = np.linalg.solve(Q, b)
fq = lambda x: 0.5 * x @ Q @ x - b @ x
x = np.zeros(n); gaps = []
for k in range(1001):
gaps.append(fq(x) - fq(xs)); g = Q @ x - b
x = x - (g @ g) / (g @ Q @ g) * g
ratio = gaps[1000] / gaps[999]
print("\nκ=800:实测每步函数值误差缩减因子 %.5f,理论上界 ((κ-1)/(κ+1))^2 = %.5f" %
(ratio, ((kappa - 1) / (kappa + 1))**2))
print("1000 步后 f-f*(初值归一化为 1): 实测 %.1e;按平方因子 %.4f;按未平方因子 %.4f" %
(gaps[1000] / gaps[0], ((kappa - 1) / (kappa + 1))**2000, ((kappa - 1) / (kappa + 1))**1000))
# ---------- Part C: 量化应用 —— 学生 t 分布极大似然:BFGS vs 最速下降,均用 Wolfe 线搜索 ----------
r = 0.01 * rng.standard_t(df=4, size=3000) + 0.0005 # 模拟厚尾日收益
def negll(th): # th = (mu, log sigma, log nu)
mu, s, eta = th; sig, nu = np.exp(s), np.exp(eta); z = (r - mu) / sig
ll = gammaln((nu + 1) / 2) - gammaln(nu / 2) - 0.5 * np.log(nu * np.pi) - s \
- (nu + 1) / 2 * np.log1p(z**2 / nu)
return -ll.mean()
def negll_g(th):
mu, s, eta = th; sig, nu = np.exp(s), np.exp(eta); z = (r - mu) / sig
d_mu = (nu + 1) * z / (sig * (nu + z**2))
d_s = -1 + (nu + 1) * z**2 / (nu + z**2)
d_nu = 0.5 * digamma((nu + 1) / 2) - 0.5 * digamma(nu / 2) - 0.5 / nu \
- 0.5 * np.log1p(z**2 / nu) + (nu + 1) * z**2 / (2 * nu * (nu + z**2))
return -np.array([d_mu.mean(), d_s.mean(), (nu * d_nu).mean()])
th0 = np.array([np.mean(r), np.log(np.std(r)), np.log(10.0)]) # 矩估计作初值
e = 1e-6; fd = np.array([(negll(th0 + e * np.eye(3)[i]) - negll(th0 - e * np.eye(3)[i])) / (2 * e) for i in range(3)])
print("\n梯度检查:解析与中心差分最大差 %.1e" % np.abs(fd - negll_g(th0)).max())
def optimize(method, th, tol=1e-6, maxit=20000):
H = np.eye(3); g = negll_g(th); unit = 0
for k in range(maxit):
if np.linalg.norm(g) < tol: return k, th, unit
p = -H @ g if method == "bfgs" else -g
a = line_search(negll, negll_g, th, p, gfk=g, c1=1e-4, c2=0.9)[0]
if a is None: a = backtracking(negll, th, negll(th), g, p)
unit += (a == 1.0)
th_new = th + a * p; g_new = negll_g(th_new)
s, y = th_new - th, g_new - g
if method == "bfgs" and y @ s > 0: # BFGS 逆更新 (2.20)
if k == 0: H = (y @ s) / (y @ y) * np.eye(3) # 首步对 H0 做尺度化(见第 8 章)
rho = 1 / (y @ s); I = np.eye(3)
H = (I - rho * np.outer(s, y)) @ H @ (I - rho * np.outer(y, s)) + rho * np.outer(s, s)
th, g = th_new, g_new
return maxit, th, unit
for m in ["bfgs", "sd"]:
k, th, unit = optimize(m, th0.copy())
print("%-4s 迭代 %5d 次,单位步长被接受 %4d 次,||g||=%.1e,-logL/T=%.6f;mu=%.5f sigma=%.5f nu=%.2f" %
(m, k, unit, np.linalg.norm(negll_g(th)), negll(th), th[0], np.exp(th[1]), np.exp(th[2])))
运行输出:
x0=(1.2, 1.2) 最速下降 18691 步 -> [1. 1.] | 牛顿 8 步 -> [1. 1.],最后 5 步步长 [1.0, 1.0, 1.0, 1.0, 1.0]
x0=(-1.2, 1.0) 最速下降 19435 步 -> [1. 1.] | 牛顿 21 步 -> [1. 1.],最后 5 步步长 [1.0, 1.0, 1.0, 1.0, 1.0]
κ=800:实测每步函数值误差缩减因子 0.99272,理论上界 ((κ-1)/(κ+1))^2 = 0.99501
1000 步后 f-f*(初值归一化为 1): 实测 4.9e-05;按平方因子 0.0067;按未平方因子 0.0821
梯度检查:解析与中心差分最大差 2.1e-10
bfgs 迭代 24 次,单位步长被接受 22 次,||g||=1.2e-08,-logL/T=-2.947158;mu=0.00049 sigma=0.00968 nu=3.87
sd 迭代 20000 次,单位步长被接受 0 次,||g||=3.1e-02,-logL/T=-2.935129;mu=0.00044 sigma=0.01105 nu=8.23
解读
- A:在 Rosenbrock 函数弯曲的山谷里,最速下降需要近 2 万步,牛顿法只要 8 步和 21 步;牛顿法最后几步的步长全是 1——这正是定理 3.5/3.7 的预言:单位步最终被接受,进入二次收敛区。
- B:在 \(\kappa=800\) 的二次函数上,渐近的每步缩减因子 0.99272 确实不超过理论上界 0.99501。实测 1000 步后误差比 0.0067 还小,是因为这个例子的特征值几何均匀分布、初期大特征值方向的误差很快被消掉,而 (3.29) 是最坏情形的界。原书"约 0.08"的数字与 (3.29) 的平方因子不一致,对应的是未平方因子,这里两个数都给出供对照。
- C:学生 \(t\) 似然的三个参数尺度差别极大(\(\mu\) 的量级是 \(10^{-4}\),对应的曲率是 \(1/\sigma^2\sim10^4\);\(\log\nu\) 的曲率是 \(O(0.1)\)),Hessian 条件数在 \(10^5\) 量级。BFGS 在 24 步内收敛到梯度 \(10^{-8}\),其中 22 步直接接受了单位步长;最速下降 2 万步后梯度仍有 \(3\times10^{-2}\),自由度估计停在 8.23(真值约 3.9),估计结果是错的,而且它不会报错。这是实务中最危险的情形:优化器"跑完了",但给出的厚尾参数严重低估了尾部风险。
实务建议:做极大似然估计时,(1) 总是先做梯度检查;(2) 优先用 BFGS/L-BFGS 或牛顿法,不要用固定步长或最速下降;(3) 结束后检查梯度范数,而不是只看 "success" 标志;(4) 参数尺度差异大时考虑重参数化或对角尺度化。
本章小结
线搜索方法的步长要满足"充分下降 + 不太短":Armijo 条件保证每步下降量与步长成比例,曲率条件(或回溯机制)防止步长过短,合起来就是 Wolfe 条件;强 Wolfe 条件进一步把步长拉近一维驻点。Zoutendijk 定理说明,只要方向不趋于与梯度垂直,Wolfe 线搜索就保证梯度趋于零。但全局收敛不等于快:最速下降按 \(\left(\frac{\kappa-1}{\kappa+1}\right)^2\) 线性收敛,病态时几乎停滞;拟牛顿法只要沿搜索方向逼近 Hessian(Dennis–Moré 条件)就超线性收敛;牛顿法局部二次收敛。两者的前提都是"总先试单位步长"。实用的线搜索用插值和保护机制实现,推荐使用 Moré–Thuente 等成熟代码。
| 概念 | 公式/要点 |
|---|---|
| Armijo 条件 | \(f(x_k+\alpha p_k)\le f_k+c_1\alpha\nabla f_k^Tp_k\),\(c_1\approx10^{-4}\) |
| 曲率条件 | \(\nabla f(x_k+\alpha p_k)^Tp_k\ge c_2\nabla f_k^Tp_k\),\(c_2=0.9\)(牛顿/拟牛顿)或 \(0.1\)(CG) |
| 强 Wolfe | \(\vert \nabla f(x_k+\alpha p_k)^Tp_k\vert \le c_2\vert \nabla f_k^Tp_k\vert \) |
| Goldstein | \(f_k+(1-c)\alpha\nabla f_k^Tp_k\le f(x_k+\alpha p_k)\le f_k+c\alpha\nabla f_k^Tp_k\),\(c<\frac12\) |
| 方向角 | \(\cos\theta_k=-\nabla f_k^Tp_k/(\Vert \nabla f_k\Vert \Vert p_k\Vert )\) |
| Zoutendijk | \(\sum_k\cos^2\theta_k\Vert \nabla f_k\Vert ^2<\infty\) |
| 条件数有界 ⇒ 方向好 | \(\Vert B_k\Vert \Vert B_k^{-1}\Vert \le M\Rightarrow\cos\theta_k\ge1/M\) |
| 最速下降收敛率 | \(\Vert x_{k+1}-x^*\Vert _Q^2\le\left(\frac{\kappa-1}{\kappa+1}\right)^2\Vert x_k-x^*\Vert _Q^2\) |
| Dennis–Moré | \(\frac{\Vert (B_k-\nabla^2f^*)p_k\Vert }{\Vert p_k\Vert }\to0\) ⇔ 超线性 |
| 牛顿法 | \(\Vert x_{k+1}-x^*\Vert \le L\Vert \nabla^2f(x^*)^{-1}\Vert \Vert x_k-x^*\Vert ^2\) |
练习
基础
- (原书 3.2)证明:若 \(0<c_2<c_1<1\),可能不存在满足 Wolfe 条件的步长。 提示:构造一个二次函数 \(\phi\),让满足充分下降的区间和满足曲率条件的区间不相交。
- (原书 3.3)证明凸二次函数沿射线的精确步长公式 (3.39)。
- (原书 3.5)证明对非奇异矩阵 \(B\) 有 \(\|Bx\|\ge\|x\|/\|B^{-1}\|\),并由此推出:若 \(B_k\) 正定且 \(\|B_k\|\|B_k^{-1}\|\le M\),则 \(\cos\theta_k\ge1/M\)。 提示:\(-\nabla f^Tp=p^TB_kp\ge\|p\|^2/\|B_k^{-1}\|\),且 \(\|\nabla f\|=\|B_kp\|\le\|B_k\|\|p\|\)。
- 推导最速下降沿 \(-\nabla f_k\) 的精确步长 (3.25),并验证 (3.27)。
- (原书 3.6)证明:若 \(x_0-x^*\) 平行于 \(Q\) 的某个特征向量,精确线搜索最速下降一步到达解。
进阶
- (原书 3.7、3.8)推导 (3.28),证明 Kantorovich 不等式,并由此得到定理 3.3。 提示:Kantorovich 不等式可通过把 \(x\) 在特征基下展开、利用凸性证明。
- (原书 3.4)证明当 \(c\le\frac12\) 时,强凸二次函数的一维精确极小点满足 Goldstein 条件。
- (原书 3.10)证明 (3.41);并证明若 \(\alpha_0\) 不满足充分下降条件,二次插值有正曲率,且 \(\alpha_1<\frac{\alpha_0}{2(1-c_1)}\)。
- 在本章量化实战 C 中,把最速下降改为先对参数做对角尺度化(用初值处 Hessian 对角元的平方根),观察迭代次数的变化。再思考:为什么 BFGS 不需要这一步? 提示:BFGS 在仿射变换下近似不变(第 2 章练习 9),它会自动"学到"尺度。
- (原书 3.9 改编)在实战 C 中把 BFGS 的线搜索参数改为 \(c_2=0.1\),统计每次迭代的函数求值次数,与 \(c_2=0.9\) 比较,说明为什么拟牛顿法默认用宽松的 \(c_2\)。
原书推荐习题:3.1 与 3.9(必做编程:回溯最速下降/牛顿与 BFGS + 强 Wolfe 在 Rosenbrock 上的对比)、3.5(条件数有界 ⇒ 方向角有界)、3.7 与 3.8(Kantorovich 不等式与最速下降收敛率)、3.2 与 3.4(Wolfe/Goldstein 参数约束的意义)、3.10 与 3.12(插值步长的性质,自己实现线搜索时有用)。
原书对照
| 本章内容 | 原书位置 | PDF 页码 |
|---|---|---|
| 3.1 线搜索框架 (3.1)(3.2) | Ch.3 引言 | PDF p.55–56 |
| 3.2.1 Wolfe 条件、引理 3.1 | §3.1 Step Length,The Wolfe Conditions | PDF p.56–61 |
| 3.2.2 Goldstein 条件 | The Goldstein Conditions | PDF p.61–62 |
| 3.2.3 回溯 Procedure 3.1 | Sufficient Decrease and Backtracking | PDF p.61–63 |
| 3.3 Zoutendijk 定理 3.2 及推论 | §3.2 Convergence of Line Search Methods | PDF p.63–66 |
| 3.4 最速下降收敛率,定理 3.3、3.4 | §3.3 Rate of Convergence,Steepest Descent | PDF p.66–69 |
| 3.5.1 Dennis–Moré,定理 3.5、3.6 | Quasi-Newton Methods | PDF p.69–71 |
| 3.5.2 牛顿法,定理 3.7 | Newton's Method | PDF p.71–73 |
| 3.6 坐标下降 | Coordinate Descent Methods,图 3.8 | PDF p.73–75 |
| 3.7 插值、初始步长、Algorithm 3.2/3.3 | §3.4 Step-Length Selection Algorithms | PDF p.75–81 |
| 注释与参考(Akaike 分析、黄金分割等) | Notes and References | PDF p.81 |
| 习题 | Exercises 3.1–3.12 | PDF p.82–83 |
页码换算:原书正文页码 = PDF 页码减 20(已用 PDF 页眉核对;第 1–2 章减 21,第 3–11 章减 20,第 12–16 章减 19,第 17 章至附录减 18)。