元信息:《Numerical Optimization》(2nd ed., Springer),Jorge Nocedal & Stephen J. Wright;本笔记负责 PDF 第 220–444 页(原书页码约 200–424),覆盖第 8 章后半至第 15 章开头。
第 8 章 拟牛顿法(Quasi-Newton Methods)(接上一块)
8.1 BFGS 方法(续):自校正性质与实现细节(PDF p.220–221)
(接上一块)上一块已推导 BFGS 逆 Hessian 更新式 (8.16)–(8.17) 与 Hessian 形式 (8.19)。记 \(s_k=x_{k+1}-x_k\),\(y_k=\nabla f_{k+1}-\nabla f_k\),\(\rho_k=1/(y_k^Ts_k)\):
坏近似能否被纠正? 若 \(y_k^Ts_k\) 很小(但为正),由 (8.16)–(8.17) 知 \(H_{k+1}\) 会变得巨大;有限精度舍入误差是否会抹掉近似矩阵的有用信息?分析与实验结论:BFGS 具有很强的自校正性质(self-correcting properties)——若 \(H_k\) 错误估计了曲率并拖慢迭代,几步之内近似会自行修正。DFP 的自校正能力弱得多,这被认为是其实际表现较差的原因。前提:必须使用合格的线搜索——Wolfe 条件保证梯度在能让模型 (8.1) 捕获恰当曲率信息的点处采样。
对偶性:DFP 与 BFGS 互为对偶,通过互换 \(s\leftrightarrow y\)、\(B\leftrightarrow H\) 可从一个公式得到另一个。
实现要点(IMPLEMENTATION):
- 线搜索应满足 Wolfe 条件 (3.6) 或强 Wolfe 条件 (3.7),且总是先试 \(\alpha_k=1\),因为在一定条件下单位步长最终总会被接受,从而获得超线性收敛。计算经验表明较不精确的线搜索更省函数求值;常用 \(c_1=10^{-4}\)、\(c_2=0.9\)。
- 初始矩阵 \(H_0=\beta I\) 没有通用的好选法。\(\beta\) 太大则首步 \(p_0=-\beta g_0\) 过长,线搜索要很多次函数求值。有些软件让用户给定首步范数 \(\delta\),取 \(H_0=\delta\|g_0\|^{-1}I\)。
- 有效的启发式:先用 \(H_0=I\) 算出第一步,但在第一次 BFGS 更新之前把 \(H_0\) 重设为
\[H_0\leftarrow\frac{y_k^Ts_k}{y_k^Ty_k}I. \quad (8.20)\]解释:设平均 Hessian \(\bar G_k=\int_0^1\nabla^2 f(x_k+\tau\alpha_kp_k)d\tau\)((8.11))正定,则 \(y_k=\bar G_ks_k\)((8.12)),存在平方根 \(\bar G_k^{1/2}\);令 \(z_k=\bar G_k^{1/2}s_k\),\[\frac{y_k^Ts_k}{y_k^Ty_k}=\frac{z_k^Tz_k}{z_k^T\bar G_kz_k}. \quad (8.21)\]其倒数是 \(\bar G_k\) 的一个 Rayleigh 商,近似 \(\nabla^2 f(x_k)\) 的某个特征值;因而 (8.21) 近似 \([\nabla^2f(x_k)]^{-1}\) 的特征值,使 \(H_0\) 的尺度与真实逆 Hessian 相当。其他缩放因子也可用,但这一个实践中最成功。
- 用 \(B_k\) 的实现:不显式存 \(B_k\),而存其 Cholesky 型分解 \(L_kD_kL_k^T\),可由 (8.19) 推出 \(O(n^2)\) 的因子更新公式;解 \(B_kp_k=-\nabla f_k\) 也只需 \(O(n^2)\)(两次三角回代 + 对角)。优点是可在 \(D_k\) 对角元过小时修正以防不稳定;但经验上无实质优势,作者偏好 Algorithm 8.1(直接更新 \(H_k\))。
- 非 Wolfe 线搜索的问题:若用 Armijo 回溯(先试 \(\alpha=1\) 再缩小直到满足充分下降 (3.6a)),不能保证曲率条件 \(y_k^Ts_k>0\)((8.7)),因为满足它可能需要大于 1 的步长。有些实现在 \(y_k^Ts_k\) 为负或接近零时跳过更新(\(H_{k+1}=H_k\)),不推荐:跳过可能过于频繁,使 \(H_k\) 无法捕获重要曲率。第 18 章的**阻尼 BFGS(damped BFGS)**是更好的对策。
8.2 SR1 方法(The SR1 Method)(PDF p.222–227)
推导。BFGS/DFP 是秩 2 修正。考虑保持对称且满足割线方程 \(B_{k+1}s_k=y_k\)((8.6))的**对称秩一(symmetric-rank-1, SR1)**更新 \(B_{k+1}=B_k+\sigma vv^T\),\(\sigma=\pm1\)。代入割线方程:
正定性:即使 \(B_k\) 正定,\(B_{k+1}\) 也可能不正定。早期只用线搜索时这被视为大缺陷;但在信赖域方法中,SR1 能产生不定近似反而是主要优点之一(能反映真实 Hessian 的不定性)。
主要缺点:分母可能为零。即使目标是凸二次函数,也可能某步不存在满足割线方程的对称秩一更新。三种情形:
- \((y_k-B_ks_k)^Ts_k\neq0\):唯一的秩一更新即 (8.24);
- \(y_k=B_ks_k\):唯一满足割线方程的"更新"为 \(B_{k+1}=B_k\);
- \(y_k\neq B_ks_k\) 且 \((y_k-B_ks_k)^Ts_k=0\):由 (8.23),不存在满足割线方程的对称秩一更新。
情形 3 表明秩一自由度不够,可能导致数值不稳定甚至崩溃——这又引回 BFGS(秩二,保证正定从而非奇异)。
仍关注 SR1 的理由:(i) 一个简单保护措施即可充分防止崩溃和数值不稳定;(ii) SR1 产生的矩阵往往是非常好的 Hessian 近似,常优于 BFGS;(iii) 约束问题的拟牛顿法、部分可分函数方法(第 18、9 章)中未必能保证 \(y_k^Ts_k>0\),BFGS 不宜,而不定近似恰能反映真实 Hessian 的不定性。
跳过规则(skipping rule):仅当
为何 SR1 可跳过而 BFGS 不宜跳过? \(s_k^T(y_k-B_ks_k)\approx0\) 需要向量以特定方式对齐,很少发生;发生时意味着 \(s_k^T\bar Gs_k\approx s_k^TB_ks_k\),即 \(B_k\) 在 \(s_k\) 方向的曲率已经正确,跳过无害。而 BFGS 的曲率条件在线搜索不施加 Wolfe 条件时(如步长不够长)很容易失败,跳过会频繁发生并损害近似质量。
Algorithm 8.2(SR1 信赖域方法):给定 \(x_0\)、\(B_0\)、信赖域半径 \(\Delta_0\)、容差 \(\epsilon>0\)、\(\eta\in(0,10^{-3})\)、\(r\in(0,1)\);\(k\leftarrow0\)。 当 \(\|\nabla f_k\|>\epsilon\):
- 求解子问题 \(\min_s \nabla f_k^Ts+\tfrac12 s^TB_ks\) s.t. \(\|s\|\le\Delta_k\)((8.27)),得 \(s_k\);
- 计算 \(y_k=\nabla f(x_k+s_k)-\nabla f_k\),实际下降 ared \(=f_k-f(x_k+s_k)\),预测下降 pred \(=-(\nabla f_k^Ts_k+\tfrac12s_k^TB_ks_k)\);
- 若 ared/pred \(>\eta\),\(x_{k+1}=x_k+s_k\),否则 \(x_{k+1}=x_k\);
- 半径更新:若 ared/pred \(>0.75\):\(\|s_k\|\le0.8\Delta_k\) 时 \(\Delta_{k+1}=\Delta_k\),否则 \(\Delta_{k+1}=2\Delta_k\);若 \(0.1\le\) ared/pred \(\le0.75\),\(\Delta_{k+1}=\Delta_k\);否则 \(\Delta_{k+1}=0.5\Delta_k\);
- 若 (8.26) 成立,用 (8.24) 计算 \(B_{k+1}\)(即使 \(x_{k+1}=x_k\) 也更新),否则 \(B_{k+1}=B_k\)。
选信赖域框架是因为无需把 Hessian 近似修改为充分正定。要点:沿失败方向也要更新 \(B_k\)——步差说明 \(B_k\) 在该方向近似差;不改善的话后续会反复生成类似方向并被拒绝,阻碍超线性收敛。
SR1 更新的性质。对二次函数,步长不影响更新,故设单位步长:\(p_k=-H_k\nabla f_k\),\(x_{k+1}=x_k+p_k\)((8.28)),于是 \(p_k=s_k\)。
定理 8.1:设 \(f(x)=b^Tx+\tfrac12x^TAx\),\(A\) 对称正定。对任意 \(x_0\) 和任意对称(不必正定)\(H_0\),只要对所有 \(k\) 有 \((s_k-H_ky_k)^Ty_k\neq0\),SR1 迭代 (8.25)(8.28) 至多 \(n\) 步收敛到极小点;若执行了 \(n\) 步且方向 \(p_i\) 线性无关,则 \(H_n=A^{-1}\)。
证明思路:归纳证明遗传性(hereditary property) \(H_ky_j=s_j\),\(j=0,\dots,k-1\)((8.29)),即割线方程沿所有历史方向成立。归纳步:由假设和 \(y_i=As_i\), \((s_k-H_ky_k)^Ty_j=s_k^Ty_j-y_k^T(H_ky_j)=s_k^Ty_j-y_k^Ts_j=s_k^TAs_j-s_k^TAs_j=0\)((8.30)),故 (8.25) 中修正项作用在 \(y_j\) 上为零,\(H_{k+1}y_j=H_ky_j=s_j\);再加上 \(H_{k+1}y_k=s_k\)。若 \(n\) 步线性无关,则 \(s_j=H_nAs_j\) 对 \(n\) 个无关向量成立,故 \(H_nA=I\),第 \(n\) 步是牛顿步,下一步即得解。若某步 \(s_k=\sum\xi_is_i\) 线性相关((8.31)),则 \(H_ky_k=\sum\xi_iH_kAs_i=\sum\xi_is_i=s_k\);又 \(s_k=-H_k\nabla f_k\),\(y_k=\nabla f_{k+1}-\nabla f_k\),得 \(H_k\nabla f_{k+1}=0\),由 \(H_k\) 非奇异得 \(\nabla f_{k+1}=0\),已到解。
意义:二次情形下 SR1 的遗传性与线搜索方式无关;BFGS 的类似结论只在精确线搜索下成立(见 8.3)。
定理 8.2(一般非线性函数):设 \(f\) 二阶连续可微,Hessian 在 \(x^*\) 邻域内有界且 Lipschitz 连续;\(\{x_k\}\) 为任意收敛到 \(x^*\) 的序列;(8.26) 对所有 \(k\) 成立;步 \(s_k\) 一致线性无关(uniformly linearly independent)。则 SR1 矩阵满足 \(\lim_{k\to\infty}\|B_k-\nabla^2f(x^*)\|=0\)。"一致线性无关"大致指步长不会趋于落在维数小于 \(n\) 的子空间中;该假设通常(但非总是)在实践中满足。注意这是 BFGS 不具备的性质(BFGS 矩阵一般不收敛到真 Hessian,只在步方向上收敛)。
8.3 Broyden 族(The Broyden Class)(PDF p.227–230)
Broyden 族:
定理 8.3:\(f\) 为强凸二次函数(Hessian \(A\)),迭代 \(p_k=-B_k^{-1}\nabla f_k\),\(x_{k+1}=x_k+p_k\)((8.34)),\(B_0\) 对称正定,用 \(\phi_k\in[0,1]\) 的 Broyden 公式更新。记 \(A^{1/2}B_k^{-1}A^{1/2}\)((8.35))的特征值 \(\lambda_1^k\le\dots\le\lambda_n^k\),则对所有 \(k\)
但最佳公式未必在受限族内:有分析与实验表明允许 \(\phi_k\) 取负(严格控制下)可能优于 BFGS。SR1 也属 Broyden 族,对应
Broyden 族的性质:
- 最后一项是秩一修正,由交错特征值定理(Theorem A.2),\(\phi_k<0\) 时特征值减小;继续减小 \(\phi_k\),矩阵先奇异后不定。\(B_{k+1}\) 奇异时
\[\phi_k^c=\frac{1}{1-\mu_k},\quad \mu_k=\frac{(y_k^TB_k^{-1}y_k)(s_k^TB_ks_k)}{(y_k^Ts_k)^2}.\quad(8.37)(8.38)\]由 Cauchy–Schwarz 不等式 \(\mu_k\ge1\),所以 \(\phi^c_k\le0\)。若 \(B_0\) 对称正定、\(s_k^Ty_k>0\) 且 \(\phi_k>\phi_k^c\),则所有 \(B_k\) 保持对称正定。
- 精确线搜索时,所有 \(\phi_k\ge\phi_k^c\) 的 Broyden 族方法产生相同的迭代序列(对一般非线性函数也成立):各方法产生的方向只差长度,精确线搜索找到同一个一维极小点。
定理 8.4(精确线搜索下的二次函数):Broyden 族方法用于强凸二次函数,\(B_0\) 对称正定,\(\alpha_k\) 为精确步长,\(\phi_k\ge\phi_k^c\)。则: (i) 至多 \(n\) 步收敛到解; (ii) 割线方程对所有历史方向成立:\(B_ks_j=y_j\),\(j=k-1,\dots,1\); (iii) 若 \(B_0=I\),迭代与共轭梯度法(第 5 章)完全相同;特别地搜索方向共轭:\(s_i^TAs_j=0\),\(i\neq j\); (iv) 若执行 \(n\) 次迭代,则 \(B_{n+1}=A\)。 推广:只要 Hessian 近似非奇异(不必正定)定理仍成立,故 \(\phi_k\) 可小于 \(\phi_k^c\) 只要不产生奇异矩阵;(iii) 推广为:若 \(B_0\neq I\),则等价于以 \(B_0\) 为预条件子的预条件共轭梯度法。这类结果主要具理论意义——实际的非精确线搜索使各方法表现差异很大——但这类分析指导了拟牛顿法的大部分发展。
8.4 收敛性分析(Convergence Analysis)(PDF p.231–239)
整体说明:BFGS 与 SR1 实践中非常稳健,但对一般非线性目标无法证明真正的全局收敛(即从任意起点、任意初始近似出发都趋于驻点),这在当时仍是未知的。分析要么假设目标凸,要么假设迭代满足某些性质。局部超线性收敛结果在合理假设下成立。本节记 \(G(x)=\nabla^2f(x)\),范数为欧氏范数。
BFGS 的全局收敛
假设 8.1:(1) \(f\) 二阶连续可微;(2) 水平集 \(\mathcal L=\{x:f(x)\le f(x_0)\}\) 凸,且存在 \(m,M>0\) 使
不能直接界住 \(B_k\) 的条件数,故引入**迹(trace)与行列式(determinant)**来估计最大、最小特征值。
定理 8.5:\(B_0\) 任意对称正定,\(x_0\) 满足假设 8.1,则 Algorithm 8.1(BFGS + Wolfe 线搜索)产生的 \(\{x_k\}\) 收敛到 \(x^*\)。
证明:定义 \(m_k=\frac{y_k^Ts_k}{s_k^Ts_k}\),\(M_k=\frac{y_k^Ty_k}{y_k^Ts_k}\)((8.42)),有 \(m_k\ge m\),\(M_k\le M\)((8.43))。对 (8.19) 取迹与行列式:
直观:\(\psi\) 同时惩罚过大特征值(迹)和过小特征值(\(-\ln\det\)),其有界性保证 \(\cos\theta_k\) 不能持续趋零——这就是"自校正"的数学体现。
推广:定理 8.5 已推广到整个受限 Broyden 族,除 DFP 外:对 \(\phi_k\in[0,1)\) 成立,\(\phi_k\to1\) 时论证失效,因为自校正性质被大大削弱。进一步可证收敛速率是线性的,且
BFGS 的超线性收敛
使用 Dennis–Moré 超线性收敛刻画 (3.32),适用于一般(非凸)非线性函数。
假设 8.2:Hessian 在 \(x^*\) 处 Lipschitz 连续:\(\|G(x)-G(x^*)\|\le L\|x-x^*\|\),\(x\) 在 \(x^*\) 附近。
变换:\(G_*=G(x^*)\),\(\tilde s_k=G_*^{1/2}s_k\),\(\tilde y_k=G_*^{-1/2}y_k\),\(\tilde B_k=G_*^{-1/2}B_kG_*^{-1/2}\);相应定义 \(\cos\tilde\theta_k\)、\(\tilde q_k\)、\(\tilde M_k=\frac{\|\tilde y_k\|^2}{\tilde y_k^T\tilde s_k}\)、\(\tilde m_k=\frac{\tilde y_k^T\tilde s_k}{\tilde s_k^T\tilde s_k}\)。BFGS 公式两边乘 \(G_*^{-1/2}\) 后形式不变:\(\tilde B_{k+1}=\tilde B_k-\frac{\tilde B_k\tilde s_k\tilde s_k^T\tilde B_k}{\tilde s_k^T\tilde B_k\tilde s_k}+\frac{\tilde y_k\tilde y_k^T}{\tilde y_k^T\tilde s_k}\),所以 (8.50) 的类比 (8.53) 成立。 由 \(y_k-G_*s_k=(\bar G_k-G_*)s_k\) 得 \(\tilde y_k-\tilde s_k=G_*^{-1/2}(\bar G_k-G_*)G_*^{-1/2}\tilde s_k\),结合假设 8.2,
定理 8.6:\(f\) 二阶连续可微,BFGS 迭代收敛到满足假设 8.2 的极小点 \(x^*\),且 (8.52) 成立,则 \(x_k\) 超线性收敛到 \(x^*\)。
证明要点:由三角不等式 \((1-\bar c\epsilon_k)\|\tilde s_k\|\le\|\tilde y_k\|\le(1+\bar c\epsilon_k)\|\tilde s_k\|\)((8.55));平方 (8.54) 得 \(\tilde m_k\ge1-\bar c\epsilon_k\)((8.56)),\(\tilde M_k\le\frac{1+\bar c\epsilon_k}{1-\bar c\epsilon_k}\)((8.57)),\(k\) 充分大时 \(\tilde M_k\le1+c\epsilon_k\)((8.58))。用 \(h(1/(1-x))\le0\) 得 \(\ln(1-x)\ge-x/(1-x)\),当 \(\bar c\epsilon_k<1/2\) 时 \(\ln\tilde m_k\ge-2\bar c\epsilon_k>-2c\epsilon_k\)((8.59))。于是
要点:超线性收敛只需 \(B_k\) 在步方向上逼近真 Hessian,而不需 \(B_k\to G_*\)。
SR1 的收敛分析
SR1 的理论不如 BFGS 完善:除二次函数结果外,没有类似定理 8.5 的全局结果或类似定理 8.6 的局部超线性结果。但对信赖域 SR1(Algorithm 8.2)有:
定理 8.7:设 \(x_k\) 由 Algorithm 8.2 产生,且 (c1) 迭代不终止,且停留在闭、有界、凸集 \(D\) 中,\(f\) 在 \(D\) 上二阶连续可微且有唯一驻点 \(x^*\);(c2) \(\nabla^2f(x^*)\) 正定,\(\nabla^2f\) 在 \(x^*\) 邻域 Lipschitz 连续;(c3) \(\{B_k\}\) 范数有界;(c4) 每步 (8.26) 成立。则 \(x_k\to x^*\),且
第 8 章 注释与参考(Notes and References)(PDF p.239–240)
综合参考:Dennis & Schnabel、Dennis & Moré、Fletcher;BFGS 的 Cholesky 因子更新见 Dennis & Schnabel。SR1 保护措施多种,(8.26) 依据 Conn–Gould–Toint 的分析更受青睐;Conn–Gould–Toint 与 Khalfan–Byrd–Schnabel 的实验表明 SR1(线搜索及信赖域版本)可与 BFGS 竞争;定理 8.7 证明见 Byrd–Khalfan–Schnabel。BFGS 全局收敛由 Powell 建立,Byrd–Nocedal–Yuan 推广到除 DFP 外的受限 Broyden 族。早期局部分析多基于有界退化原理(bounded deterioration principle):起点与初始 Hessian 近似都足够好时,迭代不会远离解,进而近似质量足以导出超线性收敛。
第 8 章习题(PDF p.240–241)
- 8.1 (a) 强凸时曲率条件 (8.7) 对任意 \(x_k,x_{k+1}\) 成立;(b) 构造一元函数 \(g(0)=-1\)、\(g(1)=-1/4\) 使 (8.7) 不成立。
- 8.2 强 Wolfe 第二条件 (3.7b) 蕴含曲率条件 (8.7)。
- 8.3 验证 (8.19) 与 (8.16) 互逆;8.4 用 Sherman–Morrison 证 (8.24) 与 (8.25) 互逆。
- 8.5 证明 (8.25) 后段落中的命题 (ii)(iii)。
- 8.6 对称正定矩阵存在对称正定平方根(用 \(A=UDU^T\))。
- 8.7 证 \(h(t)=1-t+\ln t\le0\)。8.8 证 \(\psi(B)=\sum(\lambda_i-\ln\lambda_i)>0\)。
- 8.9 证 \(\det(I+xy^T)=1+y^Tx\),\(\det(I+xy^T+uv^T)=(1+y^Tx)(1+v^Tu)-(x^Tv)(y^Tu)\),进而证明 (8.45)(题中写 (8.43) 系原书笔误)。
- 8.10 用迹的性质证 (8.44)。8.11 假设 8.1 下 \(\liminf\|\nabla f_k\|=0\) 蕴含整个序列收敛到 \(x^*\)。
第 8 章本章要点
- BFGS 有强自校正性质,前提是 Wolfe 线搜索;DFP 自校正弱。实现上先试 \(\alpha=1\),\(c_1=10^{-4},c_2=0.9\),首次更新前按 (8.20) 缩放 \(H_0\);不要在 \(y^Ts\le0\) 时简单跳过 BFGS 更新,应使用阻尼 BFGS。
- SR1 是满足割线方程的唯一对称秩一更新,可产生不定近似,适合信赖域;分母可能为零,用跳过规则 (8.26) 保护;失败步也要更新。
- SR1 对二次函数具有不依赖线搜索的遗传性,\(n\) 步内终止且 \(H_n=A^{-1}\);一般函数在一致线性无关步下 \(B_k\to\nabla^2f(x^*)\)。
- Broyden 族以 \(\phi_k\) 联结 BFGS(\(\phi=0\)) 和 DFP(\(\phi=1\));受限族使 \(A^{1/2}B_k^{-1}A^{1/2}\) 的特征值单调趋 1;\(\phi_k>\phi_k^c\) 保正定;精确线搜索时族内方法迭代相同并与 CG 等价。
- BFGS 全局收敛证明的核心是势函数 \(\psi(B)=\mathrm{tr}B-\ln\det B\);超线性收敛来自 Dennis–Moré 条件在变换坐标下的成立。SR1 在信赖域框架下有 \((n+1)\) 步超线性收敛。
第 8 章与量化交易的关联
- 参数估计与模型校准:最大似然(GARCH、随机波动率模型、期限结构模型)、隐含波动率曲面/SABR/Heston 参数校准,几乎都用 BFGS/L-BFGS(如
scipy.optimize.minimize(method='BFGS')、R 的optim)。理解 Wolfe 线搜索、\(H_0\) 缩放、\(y^Ts\le0\) 时的处理,有助于诊断"优化器卡住/不收敛"。 - BFGS 逆 Hessian 的误用:常见做法是把 BFGS 结束时的 \(H_k\) 当作 MLE 的协方差估计(标准误)。本章指出 BFGS 矩阵只在步方向上逼近真 Hessian(定理 8.6 只给出 Dennis–Moré 条件),一般不收敛到 \(\nabla^2f(x^*)\),所以用它算标准误可能严重失真,应在最优点用有限差分或解析 Hessian 重新计算。相比之下,SR1 在一致线性无关步下的近似更准确(定理 8.2)。
- 非凸问题:如带交易成本或非凸风险约束的组合优化、神经网络因子模型训练中的不定 Hessian 情形,SR1 + 信赖域比 BFGS 更合适。
第 8 章推荐习题
- 8.2(强 Wolfe ⇒ 曲率条件,理解为何 Wolfe 线搜索对 BFGS 必不可少);
- 8.4(Sherman–Morrison 推 SR1 逆公式,矩阵求逆引理的基本功);
- 8.7、8.8、8.9、8.10(全局收敛证明所需的迹、行列式、\(\psi\) 函数工具,矩阵行列式引理在其他场合也常用)。
第 9 章 大规模拟牛顿法与部分可分优化(Large-Scale Quasi-Newton and Partially Separable Optimization)(PDF p.242–269)
章引言(PDF p.242–244):第 8 章的拟牛顿近似通常是稠密的,存储和计算量随 \(n^2\) 增长,大 \(n\) 时不可承受。三条扩展路线:(1) 有限内存拟牛顿法(limited-memory quasi-Newton):只用若干长度为 \(n\) 的向量紧凑存储 Hessian 近似,稳健、廉价、易实现,但收敛不快;(2) 稀疏拟牛顿法:近似矩阵模仿真 Hessian 的稀疏模式,效果不佳,仅简述;(3) 利用多数大规模目标函数具有的**部分可分(partial separability)**结构的牛顿/拟牛顿法,通常收敛快且稳健,但需要目标函数的详细结构信息。
9.1 有限内存 BFGS(Limited-Memory BFGS, L-BFGS)(PDF p.244–249)
适用场景:Hessian 难以以合理代价计算或过于稠密。主要思想:只用最近若干次迭代的曲率信息构造 Hessian 近似,丢弃较早的(与当前 Hessian 关系较小)以节省存储。收敛速率一般只是线性,但常可接受。
回顾 BFGS:\(x_{k+1}=x_k-\alpha_kH_k\nabla f_k\)((9.1)),
L-BFGS 只隐式存储 \(m\) 对 \(\{s_i,y_i\}\),\(i=k-m,\dots,k-1\);每产生新的一对,删去最旧的一对。实践中 \(m\) 取 3 到 20 常令人满意。除更新策略与初始矩阵处理外,实现与 BFGS(Algorithm 8.1)相同,线搜索也相同。每次迭代选一个初始矩阵 \(H_k^0\)(可随迭代变化,这点不同于标准 BFGS),反复用 (9.2) 得
Algorithm 9.1(L-BFGS 双循环递归,two-loop recursion):
q ← ∇f_k
for i = k−1, k−2, ..., k−m:
α_i ← ρ_i s_iᵀ q
q ← q − α_i y_i
r ← H_k^0 q
for i = k−m, k−m+1, ..., k−1:
β ← ρ_i y_iᵀ r
r ← r + s_i (α_i − β)
return r (= H_k ∇f_k)
复杂度:不计 \(H_k^0q\),需 \(4mn\) 次乘法;\(H_k^0\) 为对角阵再加 \(n\) 次。优点:与 \(H^0_k\) 的乘法与其余计算隔离,\(H_k^0\) 可自由选取并随迭代变化;甚至可隐式地给定 Hessian 初始近似 \(B_k^0\),通过解 \(B_k^0r=q\) 得 \(r\)。
初始矩阵的有效选择:\(H_k^0=\gamma_kI\),
Algorithm 9.2(L-BFGS):选 \(x_0\)、整数 \(m>0\);\(k\leftarrow0\);重复:选 \(H_k^0\)(如 (9.6));用 Algorithm 9.1 算 \(p_k=-H_k\nabla f_k\);\(x_{k+1}=x_k+\alpha_kp_k\),\(\alpha_k\) 满足 Wolfe 条件;若 \(k>m\),丢弃 \(\{s_{k-m},y_{k-m}\}\);计算并保存 \(s_k\)、\(y_k\);\(k\leftarrow k+1\);直到收敛。 前 \(m-1\) 次迭代中,若初始矩阵相同且 \(H_k^0=H_0\) 不变,则与 BFGS 等价。理论上把 \(m\) 设得很大即可重现 BFGS,但当 \(m>n/2\) 时比直接 BFGS 更贵。保留最近 \(m\) 对的策略实践中效果好,尚无更好的一致策略;也可考虑保留使 \(S\) 矩阵条件最好的那些对(避免 \(s_i\) 近共线)。
表 9.1(数值实验):CUTE 测试集,终止准则 \(\|\nabla f_k\|\le10^{-5}\),比较 \(m=3,5,17,29\) 的函数/梯度求值次数 nfg 与 CPU 时间。例:DIXMAANL(\(n=1500\))nfg 为 146/134/120/125,时间 16.5/17.4/28.2/44.4;EIGENALS(\(n=110\))nfg 821/569/363/168,时间 21.5/15.7/16.2/12.5;FREUROTH(\(n=1000\))\(m=3,5\) 时 >999 次失败,\(m=17,29\) 时 69/38 次;TRIDIA(\(n=1000\))nfg 876/611/531/462,时间 46.6/41.4/84.6/127.1。结论:\(m\) 小时稳健性差;\(m\) 增大则函数求值次数下降,但每次迭代代价上升,最佳 CPU 时间常出现在较小的 \(m\);最优 \(m\) 依问题而定。
适用性评价:真 Hessian 不稀疏的大问题常首选 L-BFGS(计算并分解真 Hessian 的牛顿法不可行),甚至可能胜过 Newton–CG 等无 Hessian 牛顿法(第 6 章,用有限差分或自动微分算 Hessian-向量积);经验上比非线性共轭梯度法更快更稳健。若 Hessian 稠密但目标部分可分(9.4 节),利用该结构的方法在函数求值次数上远胜 L-BFGS,但按计算时间 L-BFGS 可能因单步便宜而更高效。主要弱点:常收敛慢,函数求值多;在高度病态问题(Hessian 特征值分布很宽)上低效。
与共轭梯度法的关系。有限内存方法源于改进非线性 CG 的尝试。Hestenes–Stiefel 非线性 CG((5.45))用 \(s_k=\alpha_kp_k\) 可写成
存储对比:Fletcher–Reeves \(3n\);Polak–Ribière \(4n\);Harwell VA14 \(6n\);CONMIN \(7n\);L-BFGS \(2mn+4n\)。CONMIN 是无记忆 BFGS 的扩展,因沿精心选择方向的自动重启等改进而胜过非线性 CG;VA14 是带重启的 PR 扩展。自上而下,函数求值效率与稳健性随存储增加而提高。
9.2 一般有限内存更新(General Limited-Memory Updating)(PDF p.249–253)
动机:L-BFGS 是更新逆 Hessian \(H_k\) 的线搜索法;信赖域法需要 Hessian 近似 \(B_k\);还希望有基于 SR1 的有限内存方法。用**紧凑(外积)表示(compact / outer product representation)**可高效实现所有常用拟牛顿公式及其逆,并用于约束优化中 Lagrange 函数 Hessian(或约化 Hessian)的近似(第 18 章)。只考虑每步持续刷新向量对的方法;另一种"存满就全部丢弃重新开始"的做法实践中效果较差。约定 \(B_k^{-1}=H_k\)。
BFGS 的紧凑表示——定理 9.1:\(B_0\) 对称正定,\(k\) 对 \(\{s_i,y_i\}_{i=0}^{k-1}\) 满足 \(s_i^Ty_i>0\),\(B_k\) 由 \(k\) 次 BFGS 更新 (8.19) 得到,则
有限内存版本:前 \(m\) 次迭代直接用定理 9.1,通常取 \(B_k^0=\delta_kI\),\(\delta_k=1/\gamma_k\)。\(k>m\) 时 \(S_k=[s_{k-m},\dots,s_{k-1}]\),\(Y_k=[y_{k-m},\dots,y_{k-1}]\)((9.14)),
SR1 矩阵——定理 9.2:对对称 \(B_0\) 用 \(\{s_i,y_i\}_{i=0}^{k-1}\) 做 \(k\) 次 SR1 更新,每次良定(\((y_i-B_is_i)^Ts_i\neq0\)),则
展开更新(Unrolling the Update):直观做法更贵。BFGS 写作 \(B_{k+1}=B_k-a_ka_k^T+b_kb_k^T\)((9.17)),\(a_k=B_ks_k/(s_k^TB_ks_k)^{1/2}\),\(b_k=y_k/(y_k^Ts_k)^{1/2}\)((9.18)),则 \(B_k=B_k^0+\sum_{i=k-m}^{k-1}(b_ib_i^T-a_ia_i^T)\)((9.19))。 Procedure 9.3(展开 BFGS 公式):for \(i=k-m,\dots,k-1\):\(b_i\leftarrow y_i/(y_i^Ts_i)^{1/2}\);\(a_i\leftarrow B_k^0s_i+\sum_{j=k-m}^{i-1}[(b_j^Ts_i)b_j-(a_j^Ts_i)a_j]\);\(a_i\leftarrow a_i/(s_i^Ta_i)^{1/2}\)。 \(a_i\) 依赖于被删除的最旧对,每次都得重算;\(b_i\) 与 \(b_j^Ts_i\) 可保留。设 \(B_k^0=I\),确定矩阵约需 \(\tfrac32m^2n\) 次运算,\(B v\) 需 \(4mn\)。乘积代价两者相同,但紧凑表示的更新只需 \(2mn\),远少于 \(\tfrac32m^2n\)。
9.3 稀疏拟牛顿更新(Sparse Quasi-Newton Updates)(PDF p.253–255)
思路:要求 \(B_k\) 与真 Hessian 有相同(或相似)稀疏模式,降低存储、或许更准确。设已知非零模式 \(\Omega=\{(i,j):[\nabla^2f(x)]_{ij}\neq0\text{ 对定义域内某 }x\}\),且 \((B_k)_{ij}=0,(i,j)\notin\Omega\)。定义 \(B_{k+1}\) 为二次规划的解:
9.4 部分可分函数(Partially Separable Functions)(PDF p.255–260)
可分:如 \(f(x)=f_1(x_1,x_3)+f_2(x_2,x_4,x_6)+f_3(x_5)\),每个变量只出现在一个函数中,可独立优化各部分,代价远低于 \(n\) 维优化。 部分可分:\(f\) 不可分但可写成若干**元素函数(element functions)**之和,每个元素函数沿大量线性无关方向不变。稀疏 Hessian 的函数都部分可分,许多 Hessian 不稀疏的函数也是。部分可分性带来经济的问题表示、高效的自动微分和有效的拟牛顿更新。最简单形式
简单例子:
内部变量(Internal Variables):分解方式往往不唯一,应选最细分解。例:最小曲面问题(习题 9.13),元素函数
9.5 不变子空间与部分可分性(Invariant Subspaces and Partial Separability)(PDF p.260–264)
定义 9.1(不变子空间):函数 \(f\) 的不变子空间 \(N_i\) 是 \(\mathbb R^n\) 中使 \(f(x+w)=f(x)\)(\(x,x+w\) 在定义域内)对所有 \(w\in N_i\) 成立的最大子空间。 例:\(f_i=x_{50}^2\)(\(n>50\),(9.34)),\(N_i=\{w:w_{50}=0\}\),维数 \(n-1\),紧凑表示 \(\phi_i(z)=z^2\),\(U_i=e_{50}^T\)。对比 \(f_i=(x_1+\dots+x_n)^2\)((9.35)):梯度、Hessian 完全稠密,但 \(N_i=\{w:e^Tw=0\}\)((9.36)),维数也是 \(n-1\);\(U_i=(1,\dots,1)\),\(\phi_i(z)=z^2\)。两者紧凑函数相同,只是 \(U_i\) 不同——说明稠密 Hessian 也可以有很好的部分可分结构。
定义 9.2(部分可分):\(f=\sum_{i=1}^{ne}f_i\) 且每个 \(f_i\) 有大不变子空间;即 \(f\) 可写为 (9.32) 形式,\(U_i\) 为 \(n_i\times n\) 且 \(n_i\ll n\),其零空间即 \(N_i\)。"大"的含义不精确(类似稀疏性的定义也依赖结构、规模、应用);实际中仅当所有不变子空间维数都接近 \(n\) 时才值得利用。\(U_i\) 选法不唯一(子空间基不唯一)。函数由复杂程序定义时,可用 Gay 提出的自动检测部分可分性的方法,已在 AMPL 建模语言及其自动微分中实现;第 7 章说明部分可分分解也提升自动微分效率。
稀疏性 vs 部分可分性:部分可分不要求稀疏(如 (9.35));反之稀疏 Hessian 的函数必部分可分,故部分可分更一般。 定理 9.3:二阶连续可微、Hessian 稀疏的函数都部分可分。具体地,若在开集 \(D\) 上 \(\frac{\partial^2f}{\partial x_i\partial x_j}=0\),则
群部分可分(Group Partial Separability):例如非线性最小二乘
9.6 部分可分函数的算法(Algorithms for Partially Separable Functions)(PDF p.264–267)
在牛顿法中利用:非精确(截断)牛顿法用 CG 近似解 \(\nabla^2f(x_k)p=-\nabla f(x_k)\)((9.40)),嵌入线搜索或信赖域框架(第 6 章)。CG 只需 Hessian-向量积,由 (9.33) 可通过 \(U_i\) 与 \(\nabla^2\phi_i\) 的运算得到,常比显式组装更省时省存储。最小曲面例:\(ne=n\) 个 \(2\times2\) 对称 \(\nabla^2\phi_i\),各存 3 个数;\(U_i\) 结构相同(两个 1、两个 −1,位置由 \(i\) 决定),无需显式存储。Hessian 约需 \(3n\) 存储,而标准格式存全 Hessian 下三角约需 \(5n\),节省 40%;其他应用可能节省更多。 代价:\(U_i\) 的映射需要内存访问和计算,效率损失取决于问题结构和计算机体系结构,有时显著。因此不一定追求最细分解,可以选易识别、捕捉主要结构的表示,甚至完全忽略部分可分性(如最小曲面问题可把元素函数当 4 变量函数)。 直接法:用多波前分解(multifrontal factorization,Duff–Reid,源自有限元)解 (9.40),部分组装 Hessian 并对子矩阵做稠密分解。LANCELOT 软件包允许在 CG 与多波前之间选择。结论:牛顿法中利用部分可分性的收益(按总计算时间)因问题与体系结构而异;拟牛顿法中利用该结构常使函数/梯度求值次数大幅减少。
部分可分函数的拟牛顿法:存储并更新各内部函数 Hessian 的近似 \(B_{[i]}\),组装
第 9 章 注释与参考(PDF p.267)
L-BFGS:Nocedal、Liu–Nocedal、Gilbert–Lemaréchal(后者讨论缩放参数的选择)。双循环递归依赖 BFGS 公式的特定形式 (9.2),其他 Broyden 族成员(SR1、DFP)尚无(或可能不存在)此类递归。约用一半存储的有限内存方法见 Siegel、Deuflhard 等、Gill–Leonard,是否优于 L-BFGS 未知;Buckley–Le Nir 结合拟牛顿与 CG 循环的方法常可与 L-BFGS 竞争。紧凑表示基于 Byrd–Nocedal–Schnabel [37],该文还给出非线性方程组 Broyden 矩阵的紧凑表示(第 11 章)。稀疏拟牛顿更新:Toint、Fletcher 等。部分可分概念由 Griewank–Toint 提出,系统论述见 Conn–Gould–Toint;定理 9.3 由 Griewank–Toint 证明。一般仿射变换不保持部分可分结构;部分可分拟牛顿法不具仿射不变性,但对保持可分性的变换不变,这不算缺陷。
第 9 章习题(PDF p.268–269)
- 9.1 编程实现 Algorithm 9.2,测试扩展 Rosenbrock 函数 \(f(x)=\sum_{i=1}^{n/2}[\alpha(x_{2i}-x_{2i-1}^2)^2+(1-x_{2i-1})^2]\)(\(\alpha=1,100\)),解 \(x^*=(1,\dots,1)\),起点 \((-1,\dots,-1)\),观察不同 \(m\) 的表现。
- 9.2 证 (9.7) 中 \(\hat H_{k+1}\) 奇异;9.3 精确线搜索下推导 (9.9)。
- 9.4 计算 (9.10) 的 \(B_kv\) 时如何并行,并与双循环递归的并行度比较。
- 9.5 基于 (9.16) 的有限内存 SR1:若 \(B_k^0\) 固定,如何用 \(Q_k=Y_k-B_0S_k\) 把存储减半。
- 9.6 计算 Hessian (9.27) 各元素,验证不计符号只有三种不同值。
- 9.7 把 \(f(x)=x_2x_3e^{x_1+x_3-x_4}+(x_2x_3)^2+(x_3-x_4)\) 写成 (9.32) 形式并给出各 \(U_i\)。
- 9.8 证 (9.28) 维数为 \(n-2\),(9.36) 为 \(n-1\)。
- 9.9 求两个函数的不变子空间并定义内部变量(个数应为 \(n-\dim N\))。
- 9.10 构造严格凸、但含凹元素函数的部分可分函数。
- 9.11 部分可分拟牛顿近似 (9.41)(9.45) 是否满足整体割线方程 \(Bs=y\)?
- 9.12(Griewank–Toint)仿射变换 \(t(x)=Ax+b\) 在线性部分为特定块对角时保持可分结构。
- 9.13 最小曲面问题:单位正方形上 \(q\times q\) 网格、\((q+1)^2\) 个格点,边界 \(4q\) 个点由给定函数确定,求其余高度使曲面面积最小;子方格 \(A_j\) 面积为 \(1/q^2\)(原文印为 \(q^2\)),面积 \(\iint_{A_j}\sqrt{1+z_x^2+z_y^2}\,dxdy\),用有限差分近似导数证明 \(f_j\) 具有 (9.26) 形式。
第 9 章本章要点
- L-BFGS 用最近 \(m\) 对 \((s,y)\) 隐式表示逆 Hessian,双循环递归 \(4mn\) 次乘法求 \(H_k\nabla f_k\);\(H_k^0=\gamma_kI\),\(\gamma_k=s^Ty/y^Ty\);\(m\) 取 3–20;须配合 Wolfe 线搜索。线性收敛,病态问题上慢。
- 无记忆 BFGS + 精确线搜索 = Hestenes–Stiefel CG(≈Polak–Ribière),L-BFGS 是 CG 的自然推广。
- 紧凑表示 (9.10)(9.15)(9.16) 把有限内存 BFGS/SR1 写成"基础矩阵 + 长窄矩阵 × \(2m\times2m\) 小矩阵 × 长窄矩阵",是 L-BFGS-B、有限内存 SQP 的基础;比展开公式省(\(2mn\) vs \(\tfrac32m^2n\))。
- 稀疏拟牛顿效果差;部分可分性(及群部分可分)比稀疏性更一般,内部变量 + 元素级 SR1 更新可以接近牛顿法的性能。
第 9 章与量化交易的关联
- 大规模参数估计:数千至数万参数的模型(大截面因子模型的非线性估计、带大量哑变量的面板 MLE、机器学习因子模型的 logistic/softmax 回归)通常用 L-BFGS;带上下界的参数(如波动率参数非负、权重上下限)常用 L-BFGS-B(
scipy.optimize.minimize(method='L-BFGS-B'))。\(m\) 的取舍、\(\gamma_k\) 缩放和 Wolfe 线搜索直接影响速度与稳定性;高度病态的目标(如未标准化的因子暴露)是 L-BFGS 变慢的常见原因,应先做变量缩放。 - 组合优化中的结构:组合方差 \(w^T\Sigma w\) 在因子模型 \(\Sigma=BFB^T+D\) 下是"低秩 + 对角"结构,与本章紧凑表示(对角 + 低秩外积)思想一致,可避免组装稠密 \(n\times n\) 协方差;(9.35) 类函数(如 \((\sum w_i)^2\) 预算约束罚项)说明稠密 Hessian 也可能只有一维有效结构。交易成本、按资产独立的冲击成本函数是典型可分项,多期优化目标通常是群部分可分的。
- 部分可分拟牛顿(LANCELOT 类)在量化日常工作中较少直接使用,属于了解即可的内容。
第 9 章推荐习题
- 9.1(亲手实现 L-BFGS 并观察 \(m\) 的影响,最有实践价值);
- 9.3(理解无记忆 BFGS 与 CG 的等价);
- 9.4、9.5(紧凑表示的计算与存储技巧);
- 9.7、9.9(练习识别不变子空间与内部变量)。
第 10 章 非线性最小二乘问题(Nonlinear Least-Squares Problems)(PDF p.270–295)
章引言(PDF p.271–273):最小二乘问题的目标形式为
10.1 背景(Background)(PDF p.273–279)
建模、回归与统计
例 10.1(药物浓度):在服药后时刻 \(t_j\) 抽血测浓度 \(y_j\),选模型
其他拟合准则:六次方和 \(\sum[y_j-\phi]^6\)、绝对值和 \(\sum|y_j-\phi(x;t_j)|\)((10.8))等,各有统计动机;最小二乘有坚实的统计基础。
最大似然解释:设偏差 \(\epsilon_j=y_j-\phi(x;t_j)\) 独立同分布,方差 \(\sigma^2\),密度 \(g_\sigma\)(模型准确且测量误差无系统成分时常成立)。似然
线性最小二乘(Linear Least-Squares Problems)
各 \(r_i\) 线性时 \(J\) 为常数:
- Cholesky / 正规方程法:计算 \(J^TJ\) 与 \(-J^Tr\);Cholesky 分解 \(J^TJ=\bar R^T\bar R\)((10.12),\(\bar R\) 为 \(n\times n\) 上三角,\(J\) 列满秩时存在);两次三角回代。常用且常有效,但**\(J^TJ\) 的条件数是 \(J\) 条件数的平方**;解的相对误差通常与条件数成正比,所以精度可能远差于直接处理 (10.10) 的方法;\(J\) 病态时舍入可能使对角出现小负数,分解失败。
- QR 分解法:正交变换不改变欧氏范数,\(\|Jx+r\|_2=\|Q^T(Jx+r)\|_2\)((10.13))。带列主元的 QR:
\[J\Pi=Q\begin{bmatrix}R\\0\end{bmatrix}=[Q_1\ Q_2]\begin{bmatrix}R\\0\end{bmatrix}=Q_1R,\quad(10.14)\]\(\Pi\) 为 \(n\times n\) 置换阵,\(Q\) 为 \(m\times m\) 正交,\(Q_1\) 为前 \(n\) 列,\(Q_2\) 为后 \(m-n\) 列,\(R\) 为 \(n\times n\) 上三角。于是\[\|Jx+r\|_2^2=\|R(\Pi^Tx)+Q_1^Tr\|_2^2+\|Q_2^Tr\|_2^2,\quad(10.15)\]第二项与 \(x\) 无关,令第一项为零:\(x^*=-\Pi R^{-1}Q_1^Tr\)(实际解三角系统 \(Rz=-Q_1^Tr\) 再置换 \(x^*=\Pi z\))。相对误差通常正比于 \(J\) 的条件数而非其平方,一般可靠。
- SVD 法:\(J=U\begin{bmatrix}S\\0\end{bmatrix}V^T=U_1SV^T\)((10.16)),\(U\) 为 \(m\times m\) 正交,\(V\) 为 \(n\times n\) 正交,\(S=\mathrm{diag}(\sigma_1\ge\dots\ge\sigma_n>0)\);\(J^TJ=VS^2V^T\)。同理
\[\|Jx+r\|_2^2=\|S(V^Tx)+U_1^Tr\|_2^2+\|U_2^Tr\|_2^2,\quad(10.17)\]\(x^*=-VS^{-1}U_1^Tr\),即\[x^*=-\sum_{i=1}^n\frac{u_i^Tr}{\sigma_i}v_i.\quad(10.18)\](抽取文本中此式无负号,按 (10.17) 推导应带负号,编写时请核对。)该式给出敏感性信息:\(\sigma_i\) 小时 \(x^*\) 对影响 \(u_i^Tr\) 的 \(r\) 或 \(J\) 的扰动特别敏感。\(J\) 接近秩亏(\(\sigma_n/\sigma_1\ll1\))时这一信息尤其有用,QR 一般不提供,有时值得为之付出 SVD 的额外代价。
三者比较:Cholesky 法在 \(m\gg n\) 且只能存 \(J^TJ\) 而不能存 \(J\) 时特别有用,但 \(J\) 秩亏或病态时须修改为对 \(J^TJ\) 对角元选主元。QR 中病态通常(非总是)表现为 \(R\) 右下角元素远小于其他元素,可修改策略求一个 \(J\) 略有扰动的邻近问题的解。SVD 对病态问题最稳健可靠。\(J\) 真正秩亏时部分 \(\sigma_i=0\),任意形如
10.2 非线性最小二乘算法(Algorithms for Nonlinear Least-Squares Problems)(PDF p.279–290)
Gauss–Newton 法
视为带线搜索的牛顿法的修改:去掉 Hessian 的二阶项,解
- 近似 \(\nabla^2f_k\approx J_k^TJ_k\)((10.21))省去了计算各残差 Hessian \(\nabla^2r_i\);若计算梯度 \(\nabla f_k=J_k^Tr_k\) 时已算 \(J_k\),近似几乎免费。
- 很多情形下 \(J^TJ\) 远比二阶项重要,Gauss–Newton 表现接近牛顿法。充分条件:各二阶项大小 \(|r_j(x)|\|\nabla^2r_j(x)\|\) 显著小于 \(J^TJ\) 的特征值——例如**小残差情形(small-residual case)**或 \(r_j\) 近似线性。实践中许多问题在解处残差小,常观察到快速局部收敛。
- 只要 \(J_k\) 满秩且 \(\nabla f_k\neq0\),\(p^{GN}\) 是下降方向: \((p_k^{GN})^T\nabla f_k=(p_k^{GN})^TJ_k^Tr_k=-(p_k^{GN})^TJ_k^TJ_kp_k^{GN}=-\|J_kp_k^{GN}\|_2^2\le0\),严格不等号除非 \(J_kp^{GN}=0\),后者等价于 \(J_k^Tr_k=\nabla f_k=0\)。
- (10.20) 是线性最小二乘问题
\[\min_p\|J_kp+r_k\|_2^2\quad(10.22)\]的正规方程(原书印作 \(f_k\),应为 \(r_k\)),因此可用上节 QR/SVD 求方向,不必显式形成 \(J_k^TJ_k\)。 另一动机:不是对 \(f\) 建二次模型,而是对向量函数建线性模型 \(r(x+p)\approx r(x)+J(x)p\),代入 \(\tfrac12\|r\|^2\) 后对 \(p\) 最小化。
通常沿 GN 方向做满足 Wolfe 条件 (3.6) 的线搜索,用第 3 章理论保证全局收敛。设 Jacobian 奇异值一致有下界:存在 \(\gamma>0\),
定理 10.1:各 \(r_j\) 在 \(\mathcal N\) 内 Lipschitz 连续可微,\(J\) 满足一致满秩条件 (10.23),GN 迭代步长满足 (3.6),则 \(\lim_{k\to\infty}J_k^Tr_k=0\)。 证明:\(J\) 连续,故存在 \(\beta>0\) 使 \(\|J(x)\|_2\le\beta\)(\(x\in\mathcal L\))。
局部收敛速度:取决于 \(J^TJ\) 相对二阶项的主导程度。记 \(H(x)=\sum r_j\nabla^2r_j\),类似第 3 章牛顿法分析 (3.35)–(3.37):
Levenberg–Marquardt(LM)法
把 GN 的线搜索换成信赖域,避免了 GN 在 \(J\) 秩亏或近秩亏时的弱点;仍忽略二阶项,故局部收敛性质与 GN 类似。LM 通常被视为第 4 章一般信赖域方法的鼻祖。球形信赖域子问题:
定理 10.2:Algorithm 4.1 中 \(\eta\in(0,\tfrac14)\),\(r_i\) 在水平集邻域内二阶连续可微,每步近似解满足
子问题的刻画:若 GN 解 \(\|p^{GN}\|<\Delta\),它就是 (10.26) 的解;否则存在 \(\lambda>0\) 使解 \(p^{LM}\) 满足 \(\|p\|=\Delta\) 且
LM 的实现:用第 4 章的求根算法找与给定 \(\Delta\) 近似匹配的 \(\lambda\)。易于保护:\(B=J^TJ\) 已半正定,只要 \(\lambda^{(\ell)}>0\),Cholesky 因子必存在。且不必朴素地做 Cholesky:对 (10.31) 的系数矩阵做 QR,
尺度问题与椭球信赖域:最小二乘问题常尺度很差(有的变量约 \(10^4\),有的约 \(10^{-6}\))。改用椭球信赖域:
大残差问题(Large-Residual Problems)
二阶项太大不能忽略,(10.26) 的二次模型不充分。有统计学者认为残差大说明模型不合格,应重建模型;但大残差常由人为错误造成的离群值(outliers)引起(读错仪表、地震读数归错事件),此时仍需求解,以便识别离群值并删除或降权。 GN 与 LM 在大残差情形表现通常较差:渐近收敛仅线性,慢于牛顿/拟牛顿的超线性。牛顿法(带信赖域或线搜索)在 \(\nabla^2r_j\) 易算时可选;拟牛顿法渐近也更快。但两者在早期迭代(尚未接近解时)可能不如 GN/LM,且牛顿法还需二阶导数代价。由于事先不知道残差大小,考虑混合算法(hybrid algorithms):残差小时像 GN/LM,残差大时切换到牛顿/拟牛顿步。
- Fletcher–Xu 方法:维护正定 \(B_k\)。若 GN 步使 \(f\) 下降达到某固定比例(如因子 5),则取该步并用 \(J_k^TJ_k\) 覆盖 \(B_k\);否则用 \(B_k\) 求方向并线搜索。两种情况都对 \(B_k\) 做类 BFGS 更新。零残差时最终总取 GN 步(二次收敛),非零残差时最终退化为 BFGS(超线性收敛)。Fletcher 书中表 6.1.2、6.1.3 显示在小、大、零残差问题上效果都好。
- 只近似二阶项:维护 \(S_k\approx\sum_jr_j(x_k)\nabla^2r_j(x_k)\),用 \(B_k=J_k^TJ_k+S_k\) 作信赖域或线搜索模型。Dennis–Gay–Welsch 算法(NL2SOL 软件)最著名。动机:理想地 \(S_{k+1}\approx\sum_jr_j(x_{k+1})\nabla^2r_j(x_{k+1})\);用 \((B_j)_{k+1}\) 代替各 \(\nabla^2r_j\) 并要求其沿刚走的一步模仿真值:\((B_j)_{k+1}(x_{k+1}-x_k)=\nabla r_j(x_{k+1})-\nabla r_j(x_k)\)(即 \(J\) 第 \(j\) 行之差)。由此得 \(S_{k+1}\) 的割线方程
\[S_{k+1}(x_{k+1}-x_k)=\sum_jr_j(x_{k+1})[(J_{k+1})_{j\cdot}-(J_k)_{j\cdot}]^T=J_{k+1}^Tr_{k+1}-J_k^Tr_{k+1}.\](原书此处最后一式印为 \(J_{k+1}^Tr_{k+1}-J_k^Tr_k\),按前一式应为 \(J_{k+1}^Tr_{k+1}-J_k^Tr_{k+1}=y^\sharp\),与下文 \(y^\sharp\) 定义一致。)再要求 \(S_{k+1}\) 对称、与 \(S_k\) 的差在某意义下最小,得\[S_{k+1}=S_k+\frac{(y^\sharp-S_ks)y^T+y(y^\sharp-S_ks)^T}{y^Ts}-\frac{(y^\sharp-S_ks)^Ts}{(y^Ts)^2}yy^T,\quad(10.38)\]\(s=x_{k+1}-x_k\),\(y=J_{k+1}^Tr_{k+1}-J_k^Tr_k\),\(y^\sharp=J_{k+1}^Tr_{k+1}-J_k^Tr_{k+1}\)。这是 DFP 的轻微变体(若 \(y^\sharp=y\) 则相同)。DGW 用 \(J_k^TJ_k+S_k\) 配合信赖域,另加改进:基本更新下 \(S_k\) 在趋近零残差解时不保证消失,会干扰超线性收敛,故更新前把右端 \(S_k\) 换成 \(\tau_kS_k\),\(\tau_k=\min\big(1,\frac{|s^Ty^\sharp|}{|s^TS_ks|}\big)\);当 GN 模型产生足够好的步时,从 Hessian 近似中略去 \(S_k\)。
大规模问题(Large-Scale Problems)
- \(n\) 小、\(m\) 大(如 \(m\sim10^6\)、\(n\sim100\)):\(J\) 可能太大不便存储,但可逐个计算 \(r_j\)、\(\nabla r_j\) 并累加 \(J^TJ=\sum_j\nabla r_j\nabla r_j^T\)、\(J^Tr=\sum_jr_j\nabla r_j\),直接解正规方程 (10.20)、(10.29);LM 中换 \(\lambda\) 时无需重算 \(J^TJ\),只需加 \(\lambda I\) 并重新分解。
- \(n,m\) 都大且 \(J\) 稀疏:精确解 (10.20)(10.29) 相对函数/梯度求值可能太贵。设计类似第 6 章非精确牛顿法的非精确 GN/LM:把 Hessian 换为 \(J_k^TJ_k\),半正定性简化了算法。除非 \(\nabla^2f(x^*)=J(x^*)^TJ(x^*)\),否则不能期望超线性收敛。
- Toint:把 GN 与 DGW 思想和第 9 章部分可分思想结合,维护各 \(\nabla^2r_j\) 的紧凑表示,每步决定用之或忽略(后者即 GN 步),用 CG 解步方程。
- Wright–Holt 非精确 LM:步 \(\bar p\) 满足 \(\|(J^TJ+\lambda I)\bar p+J^Tr\|\le\eta\|J^Tr\|\),\(\eta\in(0,1)\)((10.39));不用信赖域半径,回到直接操纵 \(\lambda\) 的原始策略:\(\lambda\) 取得足够大以强制"充分下降"
\[\frac{\|r(x)\|^2-\|r(x+\bar p)\|^2}{\|r(x)\|^2-\|r(x)+J(x)\bar p\|^2-\lambda^2\|\bar p\|^2}\ge\gamma_1,\quad\gamma_1\in(0,1),\quad(10.40)\](分母中 \(\lambda^2\|\bar p\|^2\) 按原文;按标准 LM 模型应为 \(\lambda\|\bar p\|^2\),编写时可核对原始文献。)且 \(\lambda\) 不比满足该式的最小值大太多,即可证全局收敛。用 Paige–Saunders 的 LSQR 算法对 (10.31) 同时求解多个 \(\lambda\):每增加一个 \(\lambda\) 只需多存两个向量和少量向量运算,无额外矩阵-向量乘积,因此代价不比解单个最小二乘问题高多少。
10.3 正交距离回归(Orthogonal Distance Regression)(PDF p.291–293)
例 10.1 假设时间 \(t_j\) 无误差,误差只来自模型不足或 \(y_j\) 测量。若忽略 \(t_j\)(自变量)中的误差,结果可能严重失真。统计学称考虑这类误差的模型为变量含误差模型(errors-in-variables models);线性模型下的优化问题称总体最小二乘(total least squares),非线性情形称正交距离回归(orthogonal distance regression, ODR)。 为 \(t_j\) 引入扰动 \(\delta_j\)、为 \(y_j\) 引入 \(\epsilon_j\):
第 10 章 注释与参考(PDF p.293–294)
大规模线性最小二乘的例子见于结构工程、数值大地测量。Lawson–Hanson 全面讨论线性最小二乘算法(含误差分析、软件清单、带界约束 \(x\ge0\) 或线性约束 \(Ax\ge b\) 的情形);Golub–Van Loan 第 5 章综述正规方程 vs QR 的适用性;Björck 的书全面综述。\(m\sim10^6\)、\(n<100\) 的问题见于 Laue 晶体学。软件:纯 GN 因缺乏稳健性不见于产品级软件,但很多算法在 GN 步有效时采用它,如 NL2SOL(Dennis–Gay–Welsch);MINPACK 有高质量 LM 实现;ODRPACK 实现正交距离回归;都允许用户提供 Jacobian 或用有限差分计算;NAG、Harwell 库也有稳健实现。
第 10 章习题(PDF p.294–295)
- 10.1 \(m\ge n\) 时 \(J\) 列满秩 ⇔ \(J^TJ\) 非奇异。
- 10.2 \(\Pi=I\) 且对角非负时 (10.12) 的 \(\bar R\) 与 (10.14) 的 \(R\) 相同。
- 10.3 \(J\) 秩亏时最小范数解对应 (10.19) 中 \(\tau_i=0\)。
- 10.4 由 \(r_j\) 及其梯度的 Lipschitz 常数 \(L\) 推出 \(J(x)\) 与 \(\nabla f(x)\) 的 Lipschitz 常数。
- 10.5 用 SVD 表示 \(\nabla f_k^Tp^{GN}\)、\(\|p^{GN}\|\)、\(\|\nabla f_k\|\) 及 \(\cos\theta_k\),说明 \(J_k\) 秩亏时为何不保证 \(\cos\theta_k>0\)。
- 10.6 用 SVD 和 \(\lambda\) 表示 (10.29) 的解及其范数平方,证明 \(\lambda\to0\) 时趋于最小范数 GN 解。
- 10.7 从 (10.46) 消去 \(p_\delta\),得到的 \(n\times n\) 系统是哪个线性最小二乘问题的正规方程。
第 10 章本章要点
- \(\nabla f=J^Tr\),\(\nabla^2f=J^TJ+\sum r_j\nabla^2r_j\);小残差或近线性时 \(J^TJ\) 主导。正态 i.i.d. 误差下最小二乘 = 最大似然。
- 线性最小二乘三法:正规方程/Cholesky(快,但条件数平方)、QR(一般首选)、SVD(最稳健,给出敏感性信息,秩亏时取最小范数解或截断小奇异值)。
- Gauss–Newton = 线性化残差后解线性最小二乘;\(J\) 一致满秩时是下降方向且全局收敛(定理 10.1);局部收敛速率由 \(\|(J^TJ)^{-1}H\|\) 决定,零残差时二次收敛。
- Levenberg–Marquardt = GN + 信赖域,等价于解 \((J^TJ+\lambda I)p=-J^Tr\);实现上用 QR + Givens 高效换 \(\lambda\),用 \(D_k\) 椭球信赖域处理尺度。
- 大残差:GN/LM 仅线性收敛,用混合法(Fletcher–Xu)或 \(J^TJ+S_k\) 的二阶项拟牛顿近似(NL2SOL)。
- 自变量也有误差时用正交距离回归,利用 Jacobian 块结构,代价与普通最小二乘相近。
第 10 章与量化交易的关联
- 因子研究与回归:横截面回归(Fama–MacBeth)、时间序列因子暴露估计本质上都是线性最小二乘。因子高度共线时,正规方程 \((X^TX)^{-1}X^Ty\) 的条件数平方效应会放大数值误差,实务中应使用 QR(
numpy.linalg.lstsq内部用 SVD,statsmodels默认用伪逆/QR);SVD 的 (10.18)(10.19) 解释了为什么共线因子的系数对数据扰动极其敏感,以及截断 SVD/主成分回归、岭回归(即 LM 的 \(J^TJ+\lambda I\) 形式)为何能稳定估计。 - 模型校准与定价:收益率曲线拟合(Nelson–Siegel–Svensson)、期权隐含波动率曲面校准(SVI、SABR、Heston 对市场价格的拟合)都是非线性最小二乘,业界标准工具就是 LM(
scipy.optimize.least_squares(method='lm'/'trf')、MINPACK、Ceres)。参数尺度差异巨大(如 Heston 的 \(\kappa\) 与 \(\rho\))时,椭球信赖域 \(D_k\) 缩放或手工标准化很关键;校准误差残差大(市场报价噪声、离群报价)时 GN/LM 收敛变慢,可考虑加权或稳健损失。 - 风险建模:GARCH 等用正态似然估计时等价于加权最小二乘结构;GN 的 \(J^TJ\) 近似对应 BHHH/外积梯度估计量,可用于标准误。
- 变量含误差:解释变量本身带估计误差(如用估计出的 beta 做第二阶段回归)会导致系数衰减偏误,正交距离回归/总体最小二乘是一种处理思路,统计上更常用工具变量或误差修正,但本章给出了计算层面的方法。
- 执行与市场微观结构:冲击成本模型(如平方根律参数)的拟合也属于非线性最小二乘。
第 10 章推荐习题
- 10.1、10.3(正规方程可解性与最小范数解,理解秩亏);
- 10.5(GN 方向在秩亏时为何失效,SVD 表示);
- 10.6(LM 解随 \(\lambda\) 变化的 SVD 表达,直接对应岭回归的收缩公式 \(\sigma_i/(\sigma_i^2+\lambda)\));
- 10.7(ODR 的块消元)。
第 11 章 非线性方程组(Nonlinear Equations)(PDF p.296–333)
章引言(PDF p.297–301):许多应用不需显式优化,而是找满足给定关系的变量值。当关系为 \(n\) 个等式、\(n\) 个变量时,即求解非线性方程组
11.1 局部算法(Local Algorithms)(PDF p.301–312)
非线性方程组的牛顿法
定理 11.1(Taylor 定理):\(r\) 在凸开集 \(D\) 内连续可微,\(x,x+p\in D\),则
Algorithm 11.1(方程组牛顿法):选 \(x_0\);对 \(k=0,1,\dots\):解牛顿方程 \(J(x_k)p_k=-r(x_k)\)((11.9));\(x_{k+1}=x_k+p_k\)。 用线性模型(而非优化中的二次模型)是因为线性模型通常有解且导出快速收敛算法。联系:无约束优化的牛顿法 (2.14) 等价于对 \(\nabla f(x)=0\) 用 Algorithm 11.1;第 18 章的 SQP 等价于对一阶最优性条件 (18.3) 用它;\(J(x_k)\) 非奇异时 (11.9) 与 Gauss–Newton (10.20) 等价。 缺点:起点远离解时行为可能混乱,\(J(x_k)\) 奇异时步甚至无定义;Jacobian 所需的一阶导数可能难以获得;\(n\) 大时精确求牛顿步可能太贵;根可能退化。退化例:\(r(x)=x^2\),唯一退化根 0,从任意 \(x_0\neq0\) 出发 \(x_k=2^{-k}x_0\),只线性收敛。
定理 11.2:\(r\) 在凸开集 \(D\) 内连续可微,\(x^*\in D\) 为非退化解,则 \(x_k\) 充分接近 \(x^*\) 时 \(x_{k+1}-x^*=o(\|x_k-x^*\|)\)((11.10),局部 Q-超线性);若 \(r\) 在 \(x^*\) 附近 Lipschitz 连续可微,则 \(x_{k+1}-x^*=O(\|x_k-x^*\|^2)\)((11.11),局部 Q-二次)。 证明:\(r(x_k)=r(x_k)-r(x^*)=J(x_k)(x_k-x^*)+o(\|x_k-x^*\|)\)((11.12))。\(J(x^*)\) 非奇异,故存在球 \(B(x^*,\delta)\)((11.13))使 \(\|J(x)^{-1}\|\le\beta^*\)((11.14))。两边乘 \(J(x_k)^{-1}\):\(-p_k=(x_k-x^*)+o(\|x_k-x^*\|)\),即 \(x_{k+1}-x^*=o(\cdot)\)((11.15))。Lipschitz 时余项 \(w(x_k,x^*)=\int_0^1[J(x_k+t(x^*-x_k))-J(x_k)](x_k-x^*)dt\) 满足 \(\|w\|=O(\|x_k-x^*\|^2)\)((11.16)(11.17)),得二次收敛。
非精确牛顿法(Inexact Newton Methods)
方向满足
定理 11.3:条件同定理 11.2,\(\{x_k\}\) 由 Framework 11.2 生成,\(x_k\) 充分接近 \(x^*\) 时:(i) 若 \(\eta\) 足够小,则 Q-线性收敛;(ii) 若 \(\eta_k\to0\),Q-超线性;(iii) 若再有 \(J\) Lipschitz 且 \(\eta_k=O(\|r_k\|)\),Q-二次。 证明:写成 \(J_kp_k=-r_k+v_k\),\(\|v_k\|\le\eta_k\|r_k\|\)((11.19));于是 \(\|p_k+J_k^{-1}r_k\|\le\beta^*\eta_k\|r_k\|\)((11.20))。\(r(x)=J(x)(x-x^*)+w(x)\),\(\rho(x)=\|w(x)\|/\|x-x^*\|\to0\)((11.21));缩小 \(\delta\) 后 \(\|r(x)\|\le2\|J(x^*)\|\|x-x^*\|+o(\cdot)\le4\|J(x^*)\|\|x-x^*\|\)((11.22))。合并得
Broyden 方法(Broyden's Method)
割线法/拟牛顿法不计算 \(J\),而构造并更新近似 \(B_k\),使其沿刚走过的步模仿真 Jacobian。步:\(p_k=-B_k^{-1}r(x_k)\),\(x_{k+1}=x_k+p_k\)((11.24))。\(s_k=x_{k+1}-x_k\),\(y_k=r(x_{k+1})-r(x_k)\),由定理 11.1,
数值例:
定理 11.5:在定理 11.2 假设下,存在 \(\epsilon,\delta>0\),若 \(\|x_0-x^*\|\le\delta\) 且 \(\|B_0-J(x^*)\|\le\epsilon\)((11.31)),则 Broyden 迭代 (11.24)(11.27) 良定且 Q-超线性收敛到 \(x^*\)(不证)。 第二个条件实际难以保证。与无约束优化不同,\(B_0\) 的选择可能对性能至关重要;一些实现建议取 \(B_0=J(x_0)\) 或其有限差分近似。即使 \(J\) 稀疏,\(B_k\) 一般稠密;\(n\) 大时可用有限内存方法:\(B_k\) 以若干长度 \(n\) 向量隐式存储,用 Sherman–Morrison–Woodbury 公式 (A.56) 解 (11.28),类似第 9 章。
张量方法(Tensor Methods)
在牛顿线性模型上加一项捕获高阶非线性行为,使对退化根(特别是 \(J(x^*)\) 秩为 \(n-1\) 或 \(n-2\) 时)收敛更快更可靠(Schnabel–Frank)。模型
11.2 实用方法(Practical Methods)(PDF p.312–324)
价值函数(Merit Functions)
单位步长的牛顿或 Broyden 法只有在离解近时才保证收敛;可能分量或 Jacobian 爆炸,也可能循环(cycling):如 \(r(x)=-x^5+x^3+4x\) 有五个非退化根,从 \(x_0=1\) 出发牛顿法在 1 和 −1 之间振荡不收敛。用线搜索和信赖域可增强稳健性,但先需价值函数(merit function)——标量函数,判断新候选点是否更接近根。最常用的是平方和
线搜索方法(Line Search Methods)
对 \(f=\tfrac12\|r\|^2\) 用第 3 章线搜索。方向须为下降方向:
定理 11.6:\(J\) 在水平集 \(\mathcal L=\{x:f(x)\le f(x_0)\}\) 的邻域 \(D\) 内 Lipschitz 连续,方向满足 (11.36),步长满足 Wolfe 条件 (3.6),则 Zoutendijk 条件成立:\(\sum_{k\ge0}\cos^2\theta_k\|J_k^Tr_k\|^2<\infty\)。 证明:令 \(\beta_R=\max(\sup_D\|r(x)\|,\sup_D\|J(x)\|)\)((11.38),原书印为 \(f(x)\),应为 \(\|r(x)\|\)),则 \(\|\nabla f(y)-\nabla f(z)\|=\|J(y)^Tr(y)-J(z)^Tr(z)\|\le\|J(y)-J(z)\|\|r(y)\|+\|J(z)\|\|r(y)-r(z)\|\le(\beta_L\beta_R+\beta_R^2)\|y-z\|\),\(\nabla f\) Lipschitz,\(f\ge0\) 有下界,定理 3.2 适用。 若保证 \(\cos\theta_k\ge\delta\)(\(\delta\in(0,1)\),\(k\) 充分大)((11.39)),则 \(J_k^Tr_k\to0\);若 \(\|J(x)^{-1}\|\) 在 \(D\) 上有界,则 \(r_k\to0\),迭代趋向一个解。
牛顿方向的 \(\cos\theta_k\):良定时,\(r_k\neq0\) 时牛顿步是下降方向:\(p_k^T\nabla f=-p_k^TJ_k^Tr_k=-\|r_k\|^2<0\),且
Algorithm 11.4(方程组的线搜索牛顿法):给定 \(\delta\in(0,1)\),\(0<c_2<c_1<\tfrac12\)(按原文;Wolfe 条件通常要求 \(c_1<c_2\),此处疑为印刷问题);选 \(x_0\);每步:若牛顿步 (11.9) 满足 (11.39) 则取之,否则用 (11.42) 并选 \(\tau_k\) 使 (11.39) 成立;若 \(\alpha=1\) 满足 Wolfe 条件 (3.6) 则 \(\alpha_k=1\),否则线搜索找满足 (3.6) 的 \(\alpha_k\);\(x_{k+1}=x_k+\alpha_kp_k\)。
定理 11.7:使用牛顿方向 (11.9) 的线搜索算法产生 \(x_k\to x^*\),\(r(x^*)=0\),\(J(x^*)\) 非奇异;\(x^*\) 邻域 \(D\) 内各 \(r_i\) 二阶可微且 \(\|\nabla^2r_i\|\) 有界;单位步满足 Wolfe 条件(\(c_1<\tfrac12\))时即接受,则 Q-二次收敛。证明略(与定理 11.10 精神类似):单位步最终总被接受,方法退化为纯牛顿法,由定理 11.2 得结论。
信赖域方法(Trust-Region Methods)
最常用的做法:对 \(f=\tfrac12\|r\|_2^2\) 用 Algorithm 4.1,模型 Hessian 取 \(B_k=J_k^TJ_k\)。全局收敛由定理 4.7、4.8 得到;精确解子问题时,在 \(J\) Lipschitz 下可证快速局部收敛。模型
狗腿法(dogleg):基于 Cauchy 点
定理 11.8:\(\eta=0\)(接受所有使价值函数下降的步),\(J\) 在水平集邻域内连续,\(\|J\|\) 在 \(\mathcal L\) 上有界,近似解满足 (11.48)(11.50),则 \(\liminf\|J_k^Tr_k\|=0\)。 定理 11.9:\(\eta\in(0,\tfrac14)\),\(J\) Lipschitz,其余同上,则 \(\lim\|J_k^Tr_k\|=0\)。
定理 11.10(局部二次收敛):Algorithm 11.5 产生的 \(x_k\) 收敛到非退化解 \(x^*\),\(J\) 在 \(x^*\) 开邻域内 Lipschitz,\(k\) 充分大时子问题精确求解,则二次收敛。意义:良好设计的算法中,为全局收敛所加的措施不干扰快速局部收敛。 证明思路:证明存在 \(K\) 使此后半径不再缩小(\(\Delta_k\ge\Delta_K\)),且最终总取纯牛顿步。(a) 无论是否受约束,\(\|p_k\|\le\|J_k^{-1}r_k\|\)((11.52))。(b) \(|1-\rho_k|\le\frac{\|r_k+J_kp_k\|^2-\|r(x_k+p_k)\|^2}{\|r_k\|^2-\|r_k+J_kp_k\|^2}\)((11.53),此处分子应理解为绝对值)。(c) 分子:\(r(x_k+p_k)=(r_k+J_kp_k)+w\),\(\|w\|\le(\beta_L/2)\|p_k\|^2\),又 \(\|r_k+J_kp_k\|\le\|r_k\|\),得分子 \(\le\epsilon(x_k)\|p_k\|^2\),\(\epsilon(x_k)=f(x_k)^{1/2}\beta_L+(\beta_L/2)^2\|p_k\|^2\)((11.54));由 \(\|p_k\|\le\beta^*\|r_k\|\to0\)((11.55))得 \(\epsilon(x_k)\to0\)((11.56))。(d) 分母:取与 \(p_k\) 等长的牛顿方向步 \(\bar p_k=-\frac{\|p_k\|}{\|J_k^{-1}r_k\|}J_k^{-1}r_k\)(可行),由最优性 \(\|r_k\|^2-\|r_k+J_kp_k\|^2\ge\|r_k\|^2-\|(1-\frac{\|p_k\|}{\|J_k^{-1}r_k\|})r_k\|^2\ge\frac{\|p_k\|}{\|J_k^{-1}r_k\|}\|r_k\|^2\ge\frac1{\beta^*}\|p_k\|\|r_k\|\)((11.57))。(e) 合并:\(|1-\rho_k|\le\frac{\beta^*\epsilon(x_k)\|p_k\|^2}{\|p_k\|\|r_k\|}\le(\beta^*)^2\epsilon(x_k)\to0\)((11.58)),故最终 \(\rho_k>\tfrac14\),半径不再缩小;而牛顿步长度 \(\le\beta^*\|r_k\|\to0\),最终小于 \(\Delta_K\),总被接受为子问题的解,由定理 11.2 得二次收敛。 可把"\(x_k\to x^*\)"弱化为"\(x^*\) 是极限点之一"(此时实际蕴含 \(x_k\to x^*\),见习题 11.11)。
11.3 延拓/同伦方法(Continuation/Homotopy Methods)(PDF p.324–330)
动机:除非 \(J(x)\) 在关心区域内非奇异(常无法保证),牛顿类方法可能收敛到价值函数的局部极小而非方程组的解。延拓法直接瞄准 \(r(x)=0\) 的解,在困难情形更可能成功。思路:构造一个解显然的"容易"方程组,逐渐变为原系统,并追踪解的移动。同伦映射(homotopy map):
实用延拓法: (1) 弧长参数化 + ODE:令 \(x,\lambda\) 为沿路径弧长 \(s\) 的函数,\((x(0),\lambda(0))=(a,0)\)。\(H(x(s),\lambda(s))=0\) 对 \(s\) 求全导:
定理 11.11(Watson):\(r\) 二阶连续可微。则对几乎所有 \(a\in\mathbb R^n\),存在从 \((a,0)\) 出发的零路径(原文写作 \((0,a)\)),沿此路径 (11.61) 满秩。若该路径在 \(\lambda\in[0,1)\) 上有界,则有聚点 \((\bar x,1)\) 满足 \(r(\bar x)=0\);若 \(J(\bar x)\) 非奇异,\((a,0)\) 到 \((\bar x,1)\) 的路径弧长有限。 即除非 \(a\) 选得不走运,路径要么发散,要么通向原系统的解。
例 11.2(延拓失败):\(r(x)=x^2-1\),两个非退化解 \(\pm1\)。取 \(a=-2\):
第 11 章 注释与参考(PDF p.330–331)
非线性微分方程和积分方程离散化是非线性方程组的重要来源;也有本身有限维的情形(如配送网络中城市间运输量)。方程 \(r_i\) 体现模型中的一致性、守恒和最优性原理。应用见 Moré、Averick 等。Broyden 法收敛分析(含定理 11.5 证明)见 Dennis–Schnabel 第 8 章、Kelley 第 6 章;有限内存 Broyden 实现见 Kelley 7.3 节。信赖域方法 (11.51) 由 Duff–Nocedal–Reid 提出。
第 11 章习题(PDF p.331–333)
- 11.1 证 \(\|ss^T/s^Ts\|=1\)。
- 11.2 \(r(x)=x^q\)(\(q>2\) 整数),牛顿法 Q-线性收敛,求收敛比。
- 11.3 验证 \(r=-x^5+x^3+4x\) 从 \(x_0=1\) 出发的循环,求根并验证非退化。
- 11.4 \(\sin(5x)-x\) 的平方和价值函数有无穷多局部极小,求通式。
- 11.5 证 \(\phi(\lambda)=\|(J^TJ+\lambda I)^{-1}J^Tr\|\) 关于 \(\lambda\) 单调减(除非 \(J^Tr=0\)),用 SVD。
- 11.6 定理 11.7 假设下 \(\nabla^2f\) Lipschitz。11.7 证定理 11.3(iii)。
- 11.8 方向只在极限意义下趋于牛顿方向(\(\|p_k+J_k^{-1}r_k\|=o(\|J_k^{-1}r_k\|)\))时仍超线性收敛。
- 11.9 精确线搜索牛顿法 \(\alpha_k\to1\)。11.10 \(JJ^Tr=0\Rightarrow J^Tr=0\)。
- 11.11* 定理 11.10 中把收敛弱化为极限点,证明 \(x^*\) 是唯一极限点(提示:证 \(\|J_{k+1}^{-1}r_{k+1}\|\le\tfrac12\|J_k^{-1}r_k\|\))。
第 11 章本章要点
- 牛顿法:解 \(J_kp=-r_k\);非退化根附近超线性(\(J\) 连续)/二次(\(J\) Lipschitz)收敛;退化根仅线性。
- 非精确牛顿:\(\|r_k+J_kp_k\|\le\eta_k\|r_k\|\),\(\eta_k\) 控制收敛速度(线性/超线性/二次),可配合 GMRES 与无矩阵 \(Jd\)。
- Broyden:满足割线方程且变化最小的秩一更新,超线性收敛,\(B_0\) 的选择比优化中重要。张量法针对退化根。
- 全局化:平方和价值函数可能有非根局部极小(只能出现在 \(J\) 奇异处);线搜索需 \(\cos\theta_k\) 有下界,牛顿方向满足 \(\cos\theta\ge1/\kappa(J)\),Powell 反例说明病态时会失败,可用 \((J^TJ+\tau I)\) 修正;信赖域(狗腿、LM)全局收敛且不破坏局部二次收敛。
- 同伦/延拓:沿 \(H(x,\lambda)=\lambda r+(1-\lambda)(x-a)\) 的零路径弧长跟踪,可越过转向点,更可靠但更贵,仍可能因路径发散而失败。
第 11 章与量化交易的关联
- 定价与校准中的求根:隐含波动率反解(一维牛顿/Brent)、从债券价格自举(bootstrap)零息曲线、多曲线框架下同时拟合多个工具的贴现/远期曲线("全局求解"方式就是 \(n\) 方程 \(n\) 未知数的非线性方程组,常用牛顿或 Broyden,Jacobian 也用于风险敏感度计算)。本章关于 \(B_0\) 选择、价值函数局部极小、Powell 反例的讨论直接对应曲线构建不收敛的排查。
- 均衡与一致性条件:风险平价组合的权重满足"各资产风险贡献相等"的非线性方程组;一般均衡/资产定价模型的不动点、结构化模型中的方程(如 Merton 模型中由股权价值和波动率反解资产价值与波动率的 2×2 方程组)都适合用本章方法。
- 同伦/延拓:参数扫描(如逐步改变相关性或约束强度来跟踪最优组合或均衡解)本质上是延拓;用于校准时可从简单模型(如 Black–Scholes 参数)逐步过渡到复杂模型,减少落入错误解的概率。
- 大规模:非精确牛顿 + GMRES 适用于 PDE 定价离散化后的非线性系统(如美式期权惩罚法、非线性 HJB 方程)。
第 11 章推荐习题
- 11.2、11.3(退化根的线性收敛与牛顿循环,体会局部方法的局限);
- 11.4(价值函数的伪局部极小);
- 11.5(LM 步长随 \(\lambda\) 单调,信赖域与 \(\lambda\) 的对应关系基础);
- 11.7(非精确牛顿二次收敛证明)。
(第 11 章习题补充,PDF p.332:11.12 把延拓失败例改为 \(r(x)=x^2-1\)、\(a=\tfrac12\),证明存在从 \((\tfrac12,0)\) 到 \((1,1)\) 的零路径(原文印为 \((1,0)\)),即换起点后延拓可成功。)
第 12 章 约束优化理论(Theory of Constrained Optimization)(PDF p.333–378)
章引言(PDF p.334–338):本书第二部分讨论带约束的极小化。一般形式
局部解与全局解:约束可能排除许多局部极小,使全局最优更易选出;但也可能使问题更难。例:\(\min\|x\|_2^2\) s.t. \(\|x\|_2^2\ge1\)——无约束时唯一解 \(x=0\);加约束后任意 \(\|x\|=1\) 都是解,\(n\ge2\) 时有无穷多局部极小。例 (12.4):\(\min(x_2+100)^2+0.01x_1^2\) s.t. \(x_2-\cos x_1\ge0\)(图 12.1),无约束唯一解 \((-100,0)\);加约束后在 \((k\pi,-1)\)(\(k=\pm1,\pm3,\pm5,\dots\))附近有大量互不连通的局部解。 定义:\(x^*\) 为局部解(local solution):\(x^*\in\Omega\) 且存在邻域 \(\mathcal N\) 使 \(f(x)\ge f(x^*)\) 对 \(x\in\mathcal N\cap\Omega\) 成立;严格(强)局部解:\(x\neq x^*\) 时严格大于;孤立局部解:\(\mathcal N\cap\Omega\) 中唯一的局部极小点。文献也常用"minimizer"代替"solution",但后者更能体现约束的作用。
光滑性:可行域边界常有棱角,但通常可用多条光滑约束描述。例:菱形区域 \(\|x\|_1=|x_1|+|x_2|\le1\)((12.5),非光滑)等价于四个线性约束 \(x_1+x_2\le1\),\(x_1-x_2\le1\),\(-x_1+x_2\le1\),\(-x_1-x_2\le1\)((12.6),图 12.2),每个约束对应多面体的一条边。非光滑无约束问题也可重写为光滑约束问题:\(f(x)=\max(x^2,x)\)((12.7),在 \(x=0,1\) 有折点,解 \(x^*=0\))等价于引入人工变量 \(t\):\(\min t\) s.t. \(t\ge x\),\(t\ge x^2\)((12.8))。当 \(f\) 是若干函数的最大值、或向量函数的 1-范数/∞-范数时常用这类技巧。任何带 \(\ge\)、\(\le\) 及非零右端的不等式都可整理为 \(c_i(x)\ge0\) 形式;一般应以直观易懂的方式陈述约束。
12.1 例子(Examples)(PDF p.338–346)
术语:可行点 \(x\) 处,不等式约束 \(i\in\mathcal I\) 若 \(c_i(x)=0\) 称活跃(active),\(c_i(x)>0\) 称非活跃(inactive)。
例 12.1(单个等式约束):\(\min x_1+x_2\) s.t. \(x_1^2+x_2^2-2=0\)((12.9))。可行集为半径 \(\sqrt2\) 的圆周,解 \(x^*=(-1,-1)^T\)。圆上其他点都能沿圆移动使 \(f\) 下降(如从 \((\sqrt2,0)\) 顺时针)。在解处约束法向 \(\nabla c_1(x^*)\) 与 \(\nabla f(x^*)\) 平行:
例 12.2(单个不等式约束):\(\min x_1+x_2\) s.t. \(2-x_1^2-x_2^2\ge0\)((12.17)),可行域为圆盘,边界上 \(\nabla c_1\) 指向内部。解仍为 \((-1,-1)\),(12.10) 在 \(\lambda_1^*=\tfrac12\) 时成立,但这里乘子的符号很重要。可行保持条件变为 \(c_1(x)+\nabla c_1(x)^Td\ge0\)((12.18))。
- 情形 I(内点,\(c_1(x)>0\)):任何足够短的 \(d\) 都满足 (12.18);\(\nabla f\neq0\) 时取 \(d=-c_1(x)\nabla f/\|\nabla f\|\) 即可同时满足 (12.13)(12.18)。唯一无此方向的情形是 \(\nabla f(x)=0\)((12.19))。
- 情形 II(边界,\(c_1(x)=0\)):条件为 \(\nabla f^Td<0\)(开半空间)与 \(\nabla c_1^Td\ge0\)(闭半空间),两者不相交当且仅当
\[\nabla f(x)=\lambda_1\nabla c_1(x),\quad\lambda_1\ge0.\quad(12.20)\]若 \(\lambda_1<0\),两梯度反向,满足条件的方向构成整个开半平面(图 12.5)。 统一表述:\[\nabla_x\mathcal L(x^*,\lambda_1^*)=0,\ \lambda_1^*\ge0,\quad(12.21)\qquad \lambda_1^*c_1(x^*)=0.\quad(12.22)\](12.22) 称互补条件(complementarity condition):乘子只有在约束活跃时才可能严格为正。情形 I 中 \(c_1>0\) 迫使 \(\lambda_1^*=0\),得 \(\nabla f=0\);情形 II 中即 (12.20)。
例 12.3(两个不等式约束):\(\min x_1+x_2\) s.t. \(2-x_1^2-x_2^2\ge0\),\(x_2\ge0\)((12.23)),可行域为半圆盘,解 \((-\sqrt2,0)^T\),两约束都活跃。一阶可行下降方向满足 \(\nabla c_i(x)^Td\ge0\)(\(i=1,2\))且 \(\nabla f(x)^Td<0\)((12.24));在解处这样的 \(d\) 不存在(满足前两条的 \(d\) 落在由 \(\nabla c_1,\nabla c_2\) 定义的象限内,而这些 \(d\) 都满足 \(\nabla f^Td\ge0\))。Lagrange 函数 \(\mathcal L=f-\lambda_1c_1-\lambda_2c_2\),条件推广为
12.2 一阶最优性条件(First-Order Optimality Conditions)(PDF p.346–350)
Lagrange 函数:
定义 12.1(LICQ,线性无关约束规范):在 \(x^*\) 处活跃约束梯度集合 \(\{\nabla c_i(x^*),i\in\mathcal A(x^*)\}\) 线性无关。(此时所有活跃约束梯度都非零。)
定理 12.1(一阶必要条件):\(x^*\) 为 (12.1) 的局部解且 LICQ 在 \(x^*\) 成立,则存在 Lagrange 乘子向量 \(\lambda^*\) 使
例 12.4:菱形区域加目标:
敏感性(Sensitivity):\(\lambda_i^*\) 刻画最优值 \(f(x^*)\) 对约束 \(c_i\) 的敏感程度——\(f\) 对该约束"推或拉"的力度。非活跃约束:微小扰动后仍不活跃,\(x^*\) 仍是局部解,\(\lambda_i^*=0\) 准确反映其无关紧要。活跃约束:把 \(c_i(x)\ge0\) 扰动为 \(c_i(x)\ge-\epsilon\|\nabla c_i(x^*)\|\),设 \(\epsilon\) 足够小使活跃集不变、乘子变化不大,则 \(-\epsilon\|\nabla c_i\|=c_i(x^*(\epsilon))-c_i(x^*)\approx(x^*(\epsilon)-x^*)^T\nabla c_i\),其他活跃约束 \(0\approx(x^*(\epsilon)-x^*)^T\nabla c_j\);由 (12.30a), \(f(x^*(\epsilon))-f(x^*)\approx\sum_j\lambda_j^*(x^*(\epsilon)-x^*)^T\nabla c_j\approx-\epsilon\|\nabla c_i\|\lambda_i^*\)。取极限:
12.3 一阶条件的推导(Derivation of the First-Order Conditions)(PDF p.350–360)
这一完整证明是理解所有约束优化算法的关键。
可行序列(feasible sequence):给定可行点 \(x^*\),序列 \(\{z_k\}\) 满足 (i) \(z_k\neq x^*\);(ii) \(z_k\to x^*\);(iii) \(k\) 充分大时 \(z_k\) 可行。所有趋于 \(x\) 的可行序列集合记为 \(T(x)\)。局部解的刻画:所有可行序列在 \(k\) 充分大时 \(f(z_k)\ge f(x)\)。 极限方向(limiting direction):\(d\) 使得沿某子序列 \(S_d\),\(\frac{z_k-x}{\|z_k-x\|}\to d\)((12.34))。\(d_k=(z_k-x)/\|z_k-x\|\) 在单位球面(紧集)上,故至少有一个极限点;可行序列可有多个极限方向;从 \(S_d\) 取元素可构造只有唯一极限方向的可行序列。
例 12.5(例 12.1 再访):在非最优点 \(x=(-\sqrt2,0)^T\) 附近,可行序列 \(z_k=(-\sqrt{2-1/k^2},-1/k)^T\)((12.35))的极限方向 \(d=(0,-1)\)(与序列在 \(x\) 处相切但方向相反)。\(f=x_1+x_2\) 沿该序列递增:\(f(z_{k+1})>f(z_k)\)(\(k\ge2\)),故 \(f(z_k)<f(x)\),\(x\) 不是解。另一序列 \(z_k=(-\sqrt{2-1/k^2},1/k)^T\) 的极限方向 \((0,1)\),\(f\) 沿它递减。两者交错组合(\(k\) 为 3 的倍数取前者,否则取后者)的序列有两个极限方向 \((0,\pm1)\)。 例 12.6(例 12.2 再访):不等式约束下可行序列多得多(图 12.10)。除上述序列外,还有从圆内沿直线趋近的 \(z_k=(-1,0)^T+(1/k)w\)((12.36)),\(w_1>0\);(原书此处以单位圆为例写 \(\|z_k\|\le1\),与 (12.17) 半径 \(\sqrt2\) 不一致,属示意性质)可行条件 \((-1+w_1/k)^2+(w_2/k)^2\le1\) 在 \(k>2w_1/(w_1^2+w_2^2)\) 时成立;也可沿曲线或随机方式趋近。
定理 12.2:若 \(x^*\) 为 (12.1) 的局部解,则 \(T(x^*)\) 中所有可行序列的任一极限方向 \(d\) 都满足
用约束规范刻画极限方向。记 \(\nabla c_i^*=\nabla c_i(x^*)\),\(A^T=[\nabla c_i^*]_{i\in\mathcal A(x^*)}\)(\(A\) 的行为活跃约束梯度),\(\nabla f^*=\nabla f(x^*)\)((12.38))。
引理 12.3:(i) 若 \(d\) 是可行序列的极限方向,则
定义 12.4:给定 \(x^*\) 与活跃集,
引入 Lagrange 乘子。LICQ 下 \(F_1\) 恰为所有可行序列极限方向的正倍数,因此定理 12.2 的条件等价于"不存在 \(d\in F_1\) 使 \(\nabla f^{*T}d<0\)"。这仍难以直接验证,下面的引理(本质上是 Farkas 引理)把它变为可检验的乘子条件。
引理 12.4:不存在 \(d\in F_1\) 使 \(d^T\nabla f^*<0\),当且仅当存在 \(\lambda\in\mathbb R^m\) 使
定理 12.1 的证明:由定理 12.2,所有极限方向满足 \(d^T\nabla f^*\ge0\);由引理 12.3(LICQ),极限方向集合恰为满足 (12.39) 的单位向量;故所有满足 (12.39) 的 \(d\) 都有 \(d^T\nabla f^*\ge0\);由引理 12.4 存在满足 (12.46) 的 \(\lambda\)。定义 \(\lambda^*_i=\lambda_i\)(\(i\in\mathcal A(x^*)\)),否则为 0((12.51)),逐条验证:(12.30a) 由 (12.46) 与定义得到;(12.30b)(12.30c) 由可行性;(12.30d):活跃不等式由 (12.46) 非负,非活跃为零;(12.30e):活跃时 \(c_i=0\),非活跃时 \(\lambda_i^*=0\)。证毕。
12.4 二阶条件(Second-Order Conditions)(PDF p.361–369)
KKT 条件成立时,沿 \(F_1\) 中任意 \(w\) 移动,一阶近似要么增加(\(w^T\nabla f>0\)),要么不变(\(w^T\nabla f=0\))。二阶导数起"决胜(tiebreaking)"作用:对 \(w^T\nabla f(x^*)=0\) 的"未决"方向,考察 Lagrange 函数的曲率。本节假设 \(f,c_i\) 二阶连续可微。 临界锥(critical cone) \(F_2(\lambda^*)\):给定满足 KKT 的 \(\lambda^*\),
定理 12.5(二阶必要条件):\(x^*\) 为局部解,LICQ 成立,\(\lambda^*\) 满足 KKT,则
充分条件:与必要条件方向相反——给出保证 \(x^*\) 为局部解的条件。与必要条件相比,不需要约束规范,且不等式改为严格。 定理 12.6(二阶充分条件):可行点 \(x^*\) 处存在 \(\lambda^*\) 满足 KKT,且
- 若 \(d\notin F_2\):存在 \(j\in\mathcal A\cap\mathcal I\) 使 \(\lambda_j^*\nabla c_j^Td>0\)((12.66)),其余 \(\lambda_i^*\nabla c_i^Td\ge0\)。则 \(\lambda_j^*c_j(z_k)=\|z_k-x^*\|\lambda_j^*\nabla c_j^Td+o(\cdot)\),\(\mathcal L(z_k,\lambda^*)\le f(z_k)-\|z_k-x^*\|\lambda_j^*\nabla c_j^Td+o(\cdot)\)((12.67));又 \(\mathcal L(z_k,\lambda^*)=f(x^*)+O(\|z_k-x^*\|^2)\),故 \(f(z_k)\ge f(x^*)+\|z_k-x^*\|\lambda_j^*\nabla c_j^Td+o(\cdot)>f(x^*)\)。
- 若 \(d\in F_2\):\(f(z_k)\ge f(x^*)+\tfrac12\|z_k-x^*\|^2d^T\nabla_{xx}\mathcal Ld+o(\cdot)>f(x^*)\)。 每个 \(z_k\) 都属于某个收敛到某极限方向的子列,故 \(k\) 充分大时 \(f(z_k)>f(x^*)\)。
例 12.7(例 12.2 第三次):\(f=x_1+x_2\),\(c_1=2-x_1^2-x_2^2\),\(\mathcal L=(x_1+x_2)-\lambda_1(2-x_1^2-x_2^2)\),KKT 在 \(x^*=(-1,-1)\)、\(\lambda_1^*=\tfrac12\) 满足;\(\nabla_{xx}\mathcal L=\mathrm{diag}(2\lambda_1^*,2\lambda_1^*)=I\) 正定,满足定理 12.6,\(x^*\) 为严格局部解(实际上是凸规划,故为全局解)。 例 12.8(非凸):
二阶条件与投影 Hessian(projected Hessians):更弱但易验证的形式。
- 最简单情形:\(\lambda^*\) 唯一(如 LICQ 成立)且严格互补成立,则 \(F_2(\lambda^*)=\mathrm{Null}[\nabla c_i(x^*)^T]_{i\in\mathcal A(x^*)}=\mathrm{Null}\,A\)。取列满秩 \(Z\) 张成该零空间,\(w=Zu\),则必要条件为 \(Z^T\nabla_{xx}\mathcal L(x^*,\lambda^*)Z\succeq0\),充分条件为 \(Z^T\nabla_{xx}\mathcal LZ\succ0\)——可通过构造此矩阵并求特征值数值验证。
- \(\lambda^*\) 唯一但严格互补不成立:\(F_2\) 是平面与半空间的交,不再是子空间。定义上下界子空间 \(\underline F_2=\{d\in F_1:\nabla c_i^Td=0,\ \forall i\in\mathcal A(x^*)\}\)(最大的含于 \(F_2\) 的子空间),\(\overline F_2=\{d\in F_1:\nabla c_i^Td=0,\ i\in\mathcal E\text{ 或 }\lambda_i^*>0\}\)(最小的包含 \(F_2\) 的子空间),\(\underline F_2\subset F_2(\lambda^*)\subset\overline F_2\)((12.70))。(原书用上下划线区分两个记号。)由 (12.55) 推出 \(\underline Z^T\nabla_{xx}\mathcal L\underline Z\succeq0\);\(\overline Z^T\nabla_{xx}\mathcal L\overline Z\succ0\) 是 (12.63) 的充分条件。这些称为双侧投影 Hessian。
- 计算 \(Z\):对 \(A^T\) 做 QR:\(A^T=Q\begin{bmatrix}R\\0\end{bmatrix}=[Q_1\ Q_2]\begin{bmatrix}R\\0\end{bmatrix}=Q_1R\)((12.71)),\(R\) 非奇异时取 \(Z=Q_2\);\(R\) 奇异(活跃约束梯度线性相关)时用带列主元的 QR(见 15.2 节)。
凸规划(Convex Programs):\(f\) 凸、可行集 \(\Omega\) 凸。所有局部解都是全局解,全局解集凸(证明类似定理 2.5,习题)。 定理 12.7:若 \(c_i\)(\(i\in\mathcal E\))线性、\(-c_i\)(\(i\in\mathcal I\))凸(即 \(c_i\) 凹),则 \(\Omega\) 凸。证明:对可行 \(x_0,x_1\) 与 \(x_\tau=(1-\tau)x_0+\tau x_1\),线性得 \(c_i(x_\tau)=0\);凸性得 \(-c_i(x_\tau)\le(1-\tau)(-c_i(x_0))+\tau(-c_i(x_1))\le0\)。 推论:线性规划是凸规划的特例。本章例 12.2、12.3、12.4 是凸规划,例 12.1 不是(可行域为圆周,非凸)。
12.5 其他约束规范(Other Constraint Qualifications)(PDF p.369–372)
约束规范的作用:保证 \(F_1\)(可行集的线性化近似)与可行序列极限方向的倍数集合一致,即线性近似捕获了真实可行集在 \(x^*\) 附近的几何特征。 失败例:\(x_2\le x_1^3\),\(x_2\ge0\)((12.72),图 12.11),\(x^*=(0,0)\) 两约束都活跃,\(F_1=\{d:-d_2\ge0,d_2\ge0\}=\{d:d_2=0\}\),含 \((-1,0)\),但它显然不是任何可行序列的极限方向;唯一可能的极限方向是 \((1,0)\)。线性近似不能反映真实几何,约束规范不成立。
引理 12.8(线性约束规范):若所有活跃约束都是线性函数,则 \(w\in F_1\) 当且仅当 \(w/\|w\|\) 是某可行序列的极限方向。证明(⇒):取 \(T>0\) 使非活跃约束沿 \(x^*+tw\)(\(t\in[0,T]\))仍为正;令 \(z_k=x^*+(T/k)w\),对活跃不等式 \(c_i(z_k)=a_i^T(z_k-x^*)=\frac Tka_i^Tw\ge0\),等式约束 \(c_i(z_k)=0\),非活跃约束由 \(T\) 的选择满足;\(\frac{z_k-x^*}{\|z_k-x^*\|}=\frac{w}{\|w\|}\)。 所以"活跃约束全为线性"是另一种约束规范,与 LICQ 互不蕴含(习题)。
定义 12.5(MFCQ,Mangasarian–Fromovitz 约束规范):存在 \(w\) 使 \(\nabla c_i(x^*)^Tw>0\)(\(i\in\mathcal A(x^*)\cap\mathcal I\),注意严格不等号),\(\nabla c_i(x^*)^Tw=0\)(\(i\in\mathcal E\)),且等式约束梯度 \(\{\nabla c_i(x^*),i\in\mathcal E\}\) 线性无关。 MFCQ 弱于 LICQ:LICQ 成立时方程组 \(\nabla c_i^Tw=1\)(活跃不等式)、\(\nabla c_i^Tw=0\)(等式)由满秩有解,即为所需 \(w\);反之可构造满足 MFCQ 而不满足 LICQ 的例子(习题 12.13)。可在 MFCQ 下证明定理 12.1;MFCQ 的良好性质:它等价于满足 KKT 的乘子集合有界(LICQ 下乘子唯一,平凡有界)。 约束规范是线性近似充分的充分条件而非必要条件:例如 \(x_2\ge-x_1^2\)、\(x_2\le x_1^2\) 在 \((0,0)\) 处上述规范都不满足,但 \(F_1=\{w:w_2=0\}\) 准确反映了可行集几何。
12.6 几何观点(A Geometric Viewpoint)(PDF p.372–375)
与代数描述 \(c_i\) 无关、只依赖可行集几何的一阶条件。问题写为 \(\min f(x)\) s.t. \(x\in\Omega\)((12.73))。用附录的切锥(tangent cone) \(T_\Omega(x^*)\)(定义 A.1)与法锥(normal cone) \(N_\Omega(x^*)\)((A.22):\(N_\Omega(x)=\{g:g^Td\le0,\ \forall d\in T_\Omega(x)\}\))。 定理 12.9:若 \(x^*\) 是 \(f\) 在 \(\Omega\) 上的局部极小点,则
第 12 章 注释与参考(PDF p.375–376)
\(N\) 的闭性(引理 12.4 所需,保证投影子问题 (12.48) 有解):\(N=\{s:s=A\lambda,\ C\lambda\ge0\}\)((12.78))。LICQ(\(A\) 列满秩)下:若 \(s_k\in N\)、\(s_k\to s^*\),则 \(\lambda_k=(A^TA)^{-1}A^Ts_k\to\lambda^*=(A^TA)^{-1}A^Ts^*\),\(C\lambda^*\ge0\),\(s^*=A\lambda^*\),故 \(s^*\in N\)。不需 LICQ 的一般证明见 Mangasarian–Schumaker(借助定理 13.2(iii))。 文献:Fletcher 第 9 章(含对偶);Bertsekas 第 3 章(强调对偶与敏感性);Mangasarian 经典著作(约束规范)。KKT 条件见 Kuhn–Tucker 1951 年论文,更早由 W. Karush 1939 年未发表硕士论文独立导出;Rockafellar 关于一般(含非光滑)问题的乘子与最优性条件的综述。
第 12 章习题(PDF p.376–378)
- 12.1 问题 (12.4) 的局部解有限还是无限多?用 KKT 论证。
- 12.2 画图说明定理 12.7 的凸性条件不必要。12.3 凸规划局部解即全局解,全局解集凸。
- 12.4 把 \(\min\|v(x)\|_\infty\)、\(\min\max_iv_i(x)\) 重写为光滑约束问题。12.5* \(\min_i f_i(x)\) 能否类似光滑重写?为什么("min of"引入非凸的析取结构)。
- 12.6 单个光滑约束可描述非光滑边界(锥 \(x_3^2\ge x_1^2+x_2^2\))。
- 12.7 证 (12.14) 在 (12.10) 不成立时满足 (12.12)(12.13)。12.8 验证 (12.35) 上 \(f\) 递增。12.9 构造趋于最大点 \((1,1)\) 的可行序列。
- 12.10 (12.72) 在原点 LICQ、MFCQ 都不成立。12.11 \(x_2\ge0\)、\(x_2\le x_1^2\) 的线性近似与约束规范。
- 12.12 活跃约束线性但 LICQ 不成立的例子。12.13 \((x_1-1)^2+(x_2-1)^2\le2\)、\((x_1-1)^2+(x_2+1)^2\le2\)、\(x_1\ge0\) 在原点 MFCQ 成立而 LICQ 不成立。12.14 验证 (12.70)。
- 12.15 求半空间 \(\{a^Tx+\alpha\ge0\}\) 中欧氏范数最小的点(投影公式)。
- 12.16(Fletcher)\(\min x_1+x_2\) s.t. \(x_1^2+x_2^2=1\) 消元时开方符号选错导致错误答案。
- 12.17 KKT + LICQ ⇒ 乘子唯一。
- 12.18 抛物线 \(y=\tfrac15(x-1)^2\) 上距 \((1,2)\) 最近的点:求 KKT 点、判断 LICQ、哪些是解;直接代入消去 \(x\) 得到的无约束问题的解不是原问题的解(消元陷阱)。
- 12.19 \(\min-2x_1+x_2\) s.t. \((1-x_1)^3-x_2\ge0\)、\(x_2+0.25x_1^2-1\ge0\),在 \(x^*=(0,1)\) 检验 LICQ 与 KKT。
- 12.20 单位圆上 \(x_1x_2\) 的极小;12.21 单位圆盘上 \(x_1x_2\) 的极大。
第 12 章本章要点
- 约束问题的局部解定义限于可行邻域;约束可能产生大量孤立局部解;非光滑问题常可通过引入辅助变量(epigraph 形式)变为光滑约束问题。
- 一阶必要条件(KKT):\(\nabla f=\sum\lambda_i\nabla c_i\),可行性,不等式乘子非负,互补松弛 \(\lambda_ic_i=0\);需约束规范(LICQ、MFCQ、线性约束等)。LICQ 下乘子唯一,MFCQ 等价于乘子集有界。
- 证明链:局部最优 ⇒ 所有可行序列极限方向非下降(定理 12.2);LICQ ⇒ 极限方向集 = 线性化锥 \(F_1\)(引理 12.3,隐函数定理);Farkas 型引理 ⇒ 存在乘子(引理 12.4,投影到闭凸锥)。
- 乘子 = 敏感性(影子价格):\(df^*/d\epsilon=-\lambda_i^*\|\nabla c_i\|\);强活跃/弱活跃之分。
- 二阶条件在临界锥 \(F_2(\lambda^*)\) 上检验 Lagrange 函数 Hessian 的曲率;实用形式是投影 Hessian \(Z^T\nabla_{xx}\mathcal LZ\);充分条件不需约束规范。
- 凸规划:等式线性、不等式 \(c_i\) 凹 ⇒ 可行域凸;KKT 即全局最优。几何形式:\(-\nabla f\in N_\Omega(x^*)\)。
第 12 章与量化交易的关联
- 组合优化的理论基础:均值-方差(Markowitz)问题 \(\min w^T\Sigma w-\gamma\mu^Tw\) s.t. \(\mathbf 1^Tw=1\)、\(w\ge0\)、行业/因子暴露约束,其 KKT 条件直接给出"所有持仓资产的边际风险调整收益相等、未持仓资产的边际收益不高于阈值"这一经典结论(互补松弛 \(\lambda_iw_i=0\))。解析求解最小方差组合、切点组合都是 KKT 的直接应用。
- 影子价格(敏感性):约束乘子告诉你放松某条约束(如个股权重上限、行业偏离、跟踪误差预算、换手率上限)每单位能带来多少目标改进,是组合经理与风控讨论约束成本的核心工具;(12.33) 的尺度说明(乘子随约束缩放改变,\(\lambda_i\|\nabla c_i\|\) 不变)在比较不同约束的"昂贵程度"时必须注意。
- 凸性判断:定理 12.7 是判断问题能否交给凸优化求解器(CVXPY 等)的基本工具:线性等式 + 凹不等式(如 \(\sigma_{\max}^2-w^T\Sigma w\ge0\))是凸的;而"持仓数不超过 K"、"最小交易单位"、"\(\|w\|=1\)"等是非凸的,例 12.1、12.8 展示了非凸约束带来的伪 KKT 点。
- 非光滑重写:最大回撤/CVaR 类目标、\(\ell_1\) 交易成本、\(\max\) 型风险指标都可用 (12.8) 的辅助变量技巧变成光滑(常为线性)约束问题,这是 Rockafellar–Uryasev CVaR 线性规划的核心手法。
- 定价:在有约束的效用最大化或无套利定价中,状态价格/风险中性测度可解释为对偶乘子。
第 12 章推荐习题
- 12.4(非光滑 → 光滑约束重写,直接用于 minimax/CVaR 建模);
- 12.13、12.10(LICQ 与 MFCQ 的区别);
- 12.15(向半空间投影的闭式解,投影类算法的基本构件);
- 12.17(乘子唯一性);
- 12.18(消元陷阱,提醒不要随意代入约束)。
第 13 章 线性规划:单纯形法(Linear Programming: The Simplex Method)(PDF p.379–410)
章引言(PDF p.380–382):Dantzig 在 20 世纪 40 年代末提出单纯形法,标志着现代优化时代的开始,使经济学家能系统高效地构建和分析大型模型,并与早期数字计算机同步发展,其计算机实现受益于数值分析而持续改进。至今线性规划与单纯形法仍是使用最广的优化工具;管理、经济、金融、工程领域的从业者长期训练于构建线性模型并用单纯形软件求解。即使实际问题是非线性的,线性规划仍有吸引力:软件成熟、保证收敛到全局最优、数据的不确定性使复杂非线性模型显得过度。内点法(第 14 章)在某些问题上更快,但单纯形法的重要性在可预见的未来不会动摇。 线性规划:线性目标 + 线性约束(等式与不等式)。可行集为多面体(polytope)——凸、连通、面为平坦多边形;目标等值线为平面(图 13.1)。解可能唯一(一个顶点),也可能是整条边、整个面,甚至整个可行集。 标准形:
- \(\min c^Tx\) s.t. \(Ax\ge b\)(无界变量):引入剩余变量(surplus) \(z\):\(Ax-z=b\),\(z\ge0\)((13.2));自由变量拆分 \(x=x^+-x^-\),\(x^+=\max(x,0)\),\(x^-=\max(-x,0)\),得 \(\min[c;-c;0]^T[x^+;x^-;z]\) s.t. \([A\ -A\ -I][x^+;x^-;z]=b\),非负。
- \(x\le u\iff x+w=u,\ w\ge0\);\(Ax\le b\iff Ax+y=b,\ y\ge0\)(松弛变量(slack))。
- \(\max c^Tx\) 改为 \(\min(-c)^Tx\)。 网络流(转运、配送)问题有特殊结构,专用单纯形算法极高效,本书不讨论(参见 Ahuja–Magnanti–Orlin)。全章设 \(m<n\);否则 \(Ax=b\) 含冗余行、不可行或只定义一个点;\(m\ge n\) 时可用 QR 或 LU 化为行满秩。
13.1 最优性与对偶(Optimality and Duality)(PDF p.383–387)
最优性条件:由第 12 章理论只需一阶 KKT 条件;凸性保证其对全局最优充分(Lagrange Hessian 为零,二阶条件无信息)。定理 12.1 需要 LICQ,但约束线性时即使相关也成立(引理 12.8)。乘子分为 \(\pi\in\mathbb R^m\)(对应 \(Ax=b\))与 \(s\in\mathbb R^n\)(对应 \(x\ge0\)),Lagrange 函数
对偶问题(dual problem):
定理 13.1(线性规划对偶定理):(i) 若 (13.1) 或 (13.7) 之一有有限最优值,则另一个也有,且最优值相等。(ii) 若其一目标无界,则另一个不可行。 证明 (i):原始有有限最优解 ⇒ 由定理 12.1 存在 \((\pi,s)\) 满足 (13.4) ⇒ 由等价性 \(\pi\) 满足 (13.8) 为对偶最优 ⇒ \(x^Ts=0\) 与 (13.10) 给出 \(c^Tx=b^T\pi\);对称论证。(ii):原始无下界 ⇒ 存在方向 \(d\):\(c^Td<0\),\(Ad=0\),\(d\ge0\);若对偶可行 \(A^T\pi\le c\),左乘 \(d^T\) 得 \(0=d^TA^T\pi\le d^Tc<0\),矛盾。
敏感性分析(sensitivity analysis):求给定最优 \(x\) 的 \((\pi,s)\) 的过程常称敏感性分析。对 \(b\) 做小扰动 \(\Delta b\),若问题非退化,扰动后的 \(s,x\) 与原来零元位置相同,互补性给出 \(x^T\Delta s=\Delta x^Ts=\dots=0\)。由对偶定理 \(c^T(x+\Delta x)=(b+\Delta b)^T(\pi+\Delta\pi)\),结合 \(c^Tx=b^T\pi\)、\(A(x+\Delta x)=b+\Delta b\)、\(A^T\Delta\pi=-\Delta s\),得 \(c^T\Delta x=\Delta b^T\pi\)。取 \(\Delta b=\epsilon e_j\):
13.2 可行集的几何(Geometry of the Feasible Set)(PDF p.387–391)
假设 \(A\) 行满秩((13.12));实践中预处理阶段会去除冗余约束、消去部分变量,并通过加松弛、剩余、人工变量保证此性质。 基本可行点(basic feasible point):可行点 \(x\),存在指标子集 \(\mathcal B(x)\subset\{1,\dots,n\}\) 满足:恰含 \(m\) 个指标;\(i\notin\mathcal B(x)\Rightarrow x_i=0\);\(m\times m\) 矩阵 \(B=[A_i]_{i\in\mathcal B(x)}\)((13.13))非奇异(\(A_i\) 为 \(A\) 的第 \(i\) 列)。单纯形法的迭代点都是基本可行点;要使该策略有意义,需 (a) 问题有基本可行点;(b) 至少一个解是基本最优点(basic optimal point)。
定理 13.2(线性规划基本定理):(i) 若 (13.1) 有可行点,则有基本可行点;(ii) 若 (13.1) 有解,则至少一个解是基本最优点;(iii) 若 (13.1) 可行且有界,则有最优解。 证明 (i):取非零分量最少(\(p\) 个,设为 \(x_1..x_p\))的可行点,\(\sum_{i\le p}A_ix_i=b\)。若 \(A_1..A_p\) 线性相关,\(A_p=\sum_{i<p}A_iz_i\)((13.14)),则 \(x(\epsilon)=x+\epsilon(z_1,\dots,z_{p-1},-1,0,\dots,0)^T\)((13.15))对任意 \(\epsilon\) 满足 \(Ax(\epsilon)=b\),且 \(|\epsilon|\) 小时前 \(p\) 个分量仍为正;存在 \(\bar\epsilon\in(0,x_p]\) 使某分量变为零,得到非零分量更少的可行点,矛盾。故列线性无关、\(p\le m\);\(p<m\) 时由 \(A\) 行满秩补充 \(m-p\) 列构成非奇异 \(B\)。(ii) 类似:对最少非零分量的最优解,若列相关,\(x^*(\epsilon)\) 对小的正负 \(\epsilon\) 都可行,最优性迫使 \(c^Tz=0\),于是可构造非零更少的最优解,矛盾。(iii) 由单纯形法有限终止得到(下一节)。 术语说明:本书"basic feasible point"即标准术语"basic feasible solution(基本可行解)","basic optimal point"即"optimal basic feasible solution";作者为与全书"solution=问题的解"一致而改用。
可行多面体的顶点:顶点是不位于集合中另外两点连线上的点(图 13.2)。 定理 13.3:(13.1) 的所有基本可行点都是可行多面体 \(\{x:Ax=b,x\ge0\}\) 的顶点,反之亦然。 证明:设 \(\mathcal B=\{1..m\}\),\(x_{m+1}=\dots=x_n=0\)((13.16))。若 \(x=\alpha y+(1-\alpha)z\)(\(\alpha\in(0,1)\),\(y,z\) 可行),则 \(y_i=z_i=0\)(\(i>m\)),\(Bx_B=By_B=Bz_B=b\),\(B\) 非奇异得 \(x=y=z\),故为顶点。反之,若顶点 \(x\) 的非零分量对应列线性相关,用 (13.15) 构造 \(x(\pm\hat\epsilon)\) 都可行,\(x\) 在其连线中点,矛盾;于是同定理 13.2 可得 \(x\) 是基本可行点。代数(基本可行点)与几何(顶点)观点一致。
定义 13.1(退化线性规划):若存在至少一个非零分量少于 \(m\) 个的基本可行点,则称该 LP 退化(degenerate)。
13.3 单纯形法(The Simplex Method)(PDF p.391–395)
方法概述:迭代点都是基本可行点(顶点);多数步从一个顶点移到相邻顶点,基指标集 \(\mathcal B\) 恰好变一个分量;多数(非全部)步使 \(c^Tx\) 下降;问题无界时,某步沿一条使目标下降且可无限前进的边。核心问题是每步决定换入/换出哪个指标,这可从 KKT 条件理解。 记 \(\mathcal N=\{1..n\}\setminus\mathcal B\)((13.17)),\(N=[A_i]_{i\in\mathcal N}\),并相应划分 \(x_B,x_N,s_B,s_N,c_B,c_N\)。由 (13.4b) \(Bx_B+Nx_N=b\),取
有限终止: \(c^Tx^+=c_B^Tx_B^++c_qx_q^+=c_B^Tx_B-c_B^TB^{-1}A_qx_q^++c_qx_q^+\)((13.23));而 \(c_B^TB^{-1}=\pi^T\),\(A_q^T\pi=c_q-s_q\),故
Procedure 13.1(单纯形法一步):给定 \(\mathcal B,\mathcal N\),\(x_B=B^{-1}b\ge0\),\(x_N=0\);
- 解 \(B^T\pi=c_B\);算 \(s_N=c_N-N^T\pi\);
- 若 \(s_N\ge0\),停止(已最优);
- 选 \(q\in\mathcal N\),\(s_q<0\) 为进基指标;
- 解 \(Bt=A_q\);若 \(t\le0\),停止(问题无界);
- 比值检验(ratio test):\(x_q^+=\min_{i:t_i>0}(x_B)_i/t_i\),达到最小的基变量指标记为 \(p\);
- 更新 \(x_B^+=x_B-tx_q^+\),\(x_N^+=(0,\dots,x_q^+,\dots,0)^T\);\(\mathcal B\) 中加入 \(q\)、移出 \(p\)。 需进一步解决三点:线性代数(维护 \(B\) 的 LU 分解以求 \(\pi,t\));从多个负 \(s_q\) 中选进基指标;处理退化步(\(x_q^+=0\),\(x\) 不变)。这些对实现效率至关重要。
13.4 单纯形法中的线性代数(Linear Algebra in the Simplex Method)(PDF p.396–400)
每步解两个系统:
13.5 其他(重要)细节(Other (Important) Details)(PDF p.400–408)
定价与进基指标选择:通常有多个负 \(s_q\)。理想是选使总步数最少的序列,但缺乏全局视角,只能用短视但实用的策略,并在搜索代价与本步下降量之间权衡。
- Dantzig 规则:计算全部 \(s_N\),选最负的 \(s_q\)。由 (13.24),它使单位 \(x_q\) 增量的目标下降最大,但不保证大的下降——可能只能把 \(x_q\) 增加一点就碰到下一顶点,甚至完全不能移动。
- 部分定价(partial pricing):每次只用 \(N\) 的一个列子矩阵算 \(s_N\) 的子向量,从中选最负者;轮换子向量使每个分量不被长期忽略。
- 多重定价(multiple pricing):对 \(s_N\) 中最负的一小批(典型 10 个)指标,不仅算价格,还算 \(t=B^{-1}A_q\) 与 \(x_q^+\),选使 \(s_qx_q^+\) 最小者;随后的迭代只在该子集上进行(相当于固定其他非基变量为零求解一个约简 LP),直到子集内价格全非负再重新定价。优点:每次只处理 \(N\) 的少数列,其余可不驻内存(内存变便宜后此优势减弱)。可组合:保留上轮最有希望的指标 + 定价 \(s_N\) 的新部分(引入"新鲜血液")。
- 最陡边(steepest edge):选沿边每单位距离使 \(c^Tx\) 下降最大的方向(Dantzig 规则是每单位 \(x_q\) 变化,二者不同:\(x_q\) 的小变化可能对应沿边的大距离)。主元步的总变化 \(x^+=x+\eta_qx_q^+\),
\[\eta_q=\begin{bmatrix}-B^{-1}A_q\\e_q\end{bmatrix}=\begin{bmatrix}-t\\e_q\end{bmatrix},\quad(13.33)\]沿 \(\eta_q\) 单位步的目标变化 \(\frac{c^T\eta_q}{\|\eta_q\|}\)((13.34)),选使之最小的 \(q\)。分子 \(c^T\eta_q=s_q\) 已知;分母用 Goldfarb–Reid 递推维护:设出基列 \(A_p\) 占 \(B\) 第一列,\(B^+=B+(A_q-Be_1)e_1^T\)((13.35));记 \(\gamma_i=\|\eta_i\|^2=\|B^{-1}A_i\|^2+1\)。Sherman–Morrison:\((B^+)^{-1}=B^{-1}-\frac{(t-e_1)e_1^TB^{-1}}{e_1^Tt}\),于是 \((B^+)^{-1}A_i=B^{-1}A_i-\frac{e_1^TB^{-1}A_i}{e_1^Tt}(t-e_1)\),进而\[\gamma_i^+=\gamma_i-2\Big(\frac{e_1^TB^{-1}A_i}{e_1^Tt}\Big)A_i^TB^{-T}t+\Big(\frac{e_1^TB^{-1}A_i}{e_1^Tt}\Big)^2\gamma_q.\quad(13.36)\]解两个系统 \(B^T\hat t=t\),\(B^Tr=e_1\)((13.37)),得\[\gamma_i^+=\gamma_i-2\Big(\frac{r^TA_i}{r^TA_q}\Big)\hat t^TA_i+\Big(\frac{r^TA_i}{r^TA_q}\Big)^2\gamma_q,\quad(13.38)\]每个 \(\gamma_i\) 的更新只需两个内积 \(r^TA_i\)、\(\hat t^TA_i\)。比 Dantzig 规则多一次 \(B^T\) 系统求解。不保证长步,但实践中非常有效;Goldfarb–Forrest 的测试显示在一些超大问题上甚至优于内点法。
启动单纯形法(两阶段法):需要初始基本可行点及基(\(B\) 非奇异,\(\bar x_B\ge0\),\(\bar x_N=0\)),找它本身与解 LP 一样难。两阶段法(Phase-I/Phase-II):
- Phase I:引入人工变量 \(z\in\mathbb R^m\),
\[\min e^Tz\quad\text{s.t.}\quad Ax+Ez=b,\ (x,z)\ge0,\quad(13.39)\]\(E\) 对角,\(E_{jj}=1\)(\(b_j\ge0\))或 \(-1\)(\(b_j<0\))。初始点 \(x=0\),\(z_j=|b_j|\)((13.40))是基本可行点,初始基矩阵即 \(E\)。\(z\) 代表对 \(Ax=b\) 的违反量,目标为违反量之和。(13.39) 最优值为零当且仅当 (13.1) 可行:若 \(e^T\tilde z=0\) 则 \(\tilde z=0\),\(\tilde x\) 可行;反之可行 \(\tilde x\) 给出 \((\tilde x,0)\),目标为 0 且目标非负,故最优。若 Phase I 终止时 \(e^Tz>0\),原问题不可行。
- Phase II:从 Phase I 的最优基出发解
\[\min c^Tx\quad\text{s.t.}\quad Ax+z=b,\ x\ge0,\ 0\ge z\ge0,\quad(13.41)\]与 (13.1) 等价(任何可行点 \(z=0\));保留 \(z\) 是因为它可能还在初始基中(取值为零),可以只保留基中的 \(z\) 分量。双侧界需对单纯形法做小修改;某个 \(z\) 分量一旦出基即从问题中删除,避免反复进出。最终解 \(z^*=0\),\(x^*\) 是 (13.1) 的基本解;若最终基仍含 \(z\) 分量,利用 \(A\) 行满秩在后处理中补入 \(x\) 指标构造 (13.1) 的最优基。
- 许多问题不需要全部 \(m\) 个人工变量:已加的松弛/剩余变量可充当人工变量。
例 13.1:\(\min3x_1+x_2+x_3\) s.t. \(2x_1+x_2+x_3\le2\),\(x_1-x_2-x_3\le-1\),\(x\ge0\)。加松弛 \(x_4,x_5\):\(2x_1+x_2+x_3+x_4=2\),\(x_1-x_2-x_3+x_5=-1\)。\(x=(0,0,0,2,0)\) 满足第一个约束但不满足第二个,只需加一个人工变量 \(z_2\):
\[\min z_2\ \text{s.t.}\ 2x_1+x_2+x_3+x_4=2,\ x_1-x_2-x_3+x_5-z_2=-1,\ (x,z_2)\ge0,\quad(13.42)\]\((x,z_2)=((0,0,0,2,0),1)\) 是基本可行点,初始基矩阵 \(B=\begin{bmatrix}1&0\\0&-1\end{bmatrix}\);\(x_4\) 充当第一个约束的人工变量。
退化步与循环(Degenerate Steps and Cycling):若某 \((x_B)_i=0\) 而 \(t_i>0\),则 \(x_q^+=0\),一步也不能移动,称退化步(degenerate step)。虽然 \(x\) 不变、目标不降,但基改变了,可能为后续下降铺路。但连续退化步可能导致循环(cycling):基变化若干次后回到之前的基,无限重复不终止。原以为罕见,但在整数规划的 LP 松弛中越来越常见,实用代码必须有防循环策略。 扰动策略:把 \(b\) 扰动为 \(b(\epsilon)=b+E(\epsilon,\epsilon^2,\dots,\epsilon^m)^T\)(\(E\) 非奇异,\(\epsilon>0\) 小),则
13.6 单纯形法的定位(Where Does the Simplex Method Fit?)(PDF p.408–409)
含不等式约束(含界约束)的优化问题的根本任务是把不等式划分为在解处活跃与不活跃两类。单纯形法属于有效集方法(active set methods):显式维护活跃/非活跃指标集估计并每步更新(LP 中 \(\mathcal B\) 是"可能不活跃"指标——\(x_i\ge0\) 不活跃,\(\mathcal N\) 是"可能活跃"指标),每步只在两集合间交换一个指标。二次规划、界约束优化、非线性规划的有效集算法沿用同一策略:定义"可能活跃"集,朝把这些约束当等式的约简问题的解前进。非线性情形下许多使单纯形高效的特征不再成立:解处不再一定有至少 \(n-m\) 个界活跃;专用线性代数不再适用;约简问题未必能每步精确求解。但单纯形法仍是有效集方法的鼻祖。 复杂度:几乎所有实际问题上很高效(一般至多 \(2m\) 到 \(3m\) 次迭代),但存在病态例:Klee–Minty 构造了可行多面体有 \(2^n\) 个顶点的 \(n\) 维问题,单纯形法访问每个顶点才到最优,证明单纯形法复杂度是指数的。多年来人们寻找多项式算法:70 年代末 Khachiyan 的椭球法是多项式的,但实践中太慢;80 年代中 Karmarkar 提出从可行多面体内部逼近解的多项式算法,开启了内点法的研究热潮(第 14 章)。
第 13 章 注释与参考(PDF p.409)
Wolfe 提出另一种 Phase I:不引入人工变量,从任意满足 \(Ax=b\)、至多 \(m\) 个非零分量的 \(x\) 出发(\(x_B\) 不必全正),求解 \(\min\sum_{x_i<0}(-x_i)\) s.t. \(Ax=b\),目标为零时终止。目标只是分段线性,但仍可用单纯形法:每步重新定义成本向量 \(f_i=-1\)(\(x_i<0\))否则 0。
第 13 章习题(PDF p.410)
- 13.1 超定系统 \(Ax=b\) 的完全主元 Gauss 消元 \(PAQ=L\begin{bmatrix}U_{11}&U_{12}\\0&0\end{bmatrix}\):(a) 可行当且仅当 \(L^{-1}Pb\) 后 \(m-\bar m\) 分量为零;(b) \(\bar m=n\) 时求唯一解;(c) 前 \(\bar m\) 行构成的约简系统等价于原系统。(对应预处理中去除冗余行。)
- 13.2 把 \(\max c^Tx+d^Ty\) s.t. \(A_1x=b_1\)、\(A_2x+B_2y\le b_2\)、\(l\le y\le u\)(\(x\) 无界)化为标准形。
- 13.3 证 \(\min c^Tx\) s.t. \(Ax\ge b,x\ge0\) 的对偶是 \(\max b^T\pi\) s.t. \(A^T\pi\le c,\pi\ge0\)。
- 13.4 \(A\) 行相关时 \(B\) 奇异,无基本可行点。
- 13.5 验证 Goldfarb–Reid 公式 (13.36)。13.6 由 \(L_1U_1\) 最后一行求 \(l_{52},l_{53},l_{54},\hat w_2\)。13.7 如何用 (13.31) 高效解 \(B^+\) 系统。
第 13 章本章要点
- 任何 LP 都可化为标准形 \(\min c^Tx\), \(Ax=b\), \(x\ge0\);KKT 条件 = 原始可行 + 对偶可行 + 互补松弛,且对 LP 充分必要。
- 对偶:\(\max b^T\pi\), \(A^T\pi\le c\);弱对偶 \(c^Tx\ge b^T\pi\);强对偶(有限最优值相等);一方无界则另一方不可行。对偶变量 = 影子价格 \(\partial(c^Tx^*)/\partial b_j=\pi_j\)。
- 基本可行点 = 顶点;有解则有基本最优解,因此只需在顶点中搜索。
- 单纯形一步:解 \(B^T\pi=c_B\) 定价得检验数 \(s_N\);选负检验数进基;解 \(Bt=A_q\);比值检验选出基;目标下降 \(-s_qx_q^+\)。非退化时有限终止。
- 实现:LU 分解 + Forrest–Tomlin/Bartels–Golub 更新 + 定期重分解;定价规则(Dantzig、部分/多重定价、最陡边);两阶段法求初始基;扰动/字典序防循环。
- 单纯形法是有效集方法的原型;实际高效但最坏指数复杂度(Klee–Minty),催生多项式时间的椭球法与内点法。
第 13 章与量化交易的关联
- 组合构建中的 LP:最小化 CVaR(Rockafellar–Uryasev 形式)、\(\ell_1\)/最小绝对偏差(MAD)风险模型、带线性交易成本和换手约束的组合再平衡、指数跟踪中的 \(\ell_1\) 跟踪误差最小化,都可写成 LP。将 \(|w_i-w_i^0|\) 拆分为买入/卖出两个非负变量正是本章"\(x=x^+-x^-\)"技巧。
- 影子价格:LP 对偶变量直接给出每条约束(如行业中性、因子暴露上限、流动性上限)的边际成本,用于约束取舍与报告;敏感性分析 (13.11) 说明在非退化条件下影子价格在局部是线性的,但基改变时会跳变——解释了为何组合对约束参数的响应呈分段线性。
- 基本最优解的稀疏性:LP 最优解至多有 \(m\) 个非零分量(定理 13.2),这解释了为何纯 LP 形式的组合优化(如 CVaR 最小化)往往给出集中、稀疏的持仓。
- 无套利与定价:资产定价基本定理的离散版本就是 LP 对偶(Farkas 引理):无套利 ⇔ 存在正的状态价格向量;超额复制(super-replication)价格是一个 LP,其对偶是在风险中性测度集合上求期望最大值。
- 执行与调度:订单拆分、跨交易所路由的简化模型可用 LP;混合整数规划(如最小交易手数、持仓数量限制)的分支定界依赖 LP 松弛的反复求解,正是循环问题变得常见的场景。
- 实务中直接调用 HiGHS、Gurobi、CPLEX、
scipy.optimize.linprog(method='highs'),本章帮助理解输出(基状态、对偶值、检验数)与求解失败(不可行、无界、退化)的含义。
第 13 章推荐习题
- 13.2、13.3(标准形转换与对偶推导,建模基本功);
- 13.4(行满秩假设的作用);
- 13.5(最陡边递推,熟悉 Sherman–Morrison 在基更新中的使用);
- 13.1(预处理中冗余约束的检测)。
第 14 章 线性规划:内点法(Linear Programming: Interior-Point Methods)(PDF p.411–436)
章引言(PDF p.412–413):20 世纪 80 年代发现,把许多大型线性规划表述为非线性问题、用牛顿法等非线性算法的变体求解非常高效。这类方法要求所有迭代点严格满足不等式约束,故称内点法(interior-point methods)。到 90 年代初,**原始–对偶方法(primal–dual methods)**脱颖而出,成为最有效的实用方法,在大型问题上与单纯形法强有力竞争。 动机:单纯形法最坏情况下是指数时间。Khachiyan 椭球法在最坏情况下是多项式的,但在所有问题上都接近最坏界,不具竞争力。Karmarkar 1984 年的投影算法既多项式又有良好实际表现(最初的出色性能宣传并未完全兑现),引发大量研究:仿射尺度(affine-scaling)、对数障碍(logarithmic barrier)、势函数下降(potential-reduction)、路径跟踪(path-following)、原始–对偶、不可行内点等,均与 Karmarkar 算法及 17.2 节对数障碍法相关。 与单纯形法对比:内点法每次迭代代价高但进展大;单纯形法迭代多而便宜。单纯形沿可行多面体边界逐个检查顶点;内点法只在极限时接近边界,可从内部或外部逼近解,但从不位于边界上。 本章内容:原始–对偶内点法的基本思想(与牛顿法、同伦法的关系,中心路径,中心邻域)、路径跟踪方法、Mehrotra 预测–校正算法(当前软件的基础)。
14.1 原始–对偶方法(Primal–Dual Methods)(PDF p.413–421)
框架:标准形原始问题
中心路径(The Central Path):严格可行点构成的弧 \(\mathcal C\),以 \(\tau>0\) 参数化,每点 \((x_\tau,\lambda_\tau,s_\tau)\) 满足
Framework 14.1(原始–对偶框架):给定 \((x^0,\lambda^0,s^0)\in\mathcal F^o\);对 \(k=0,1,\dots\):解
不可行内点法(infeasible-interior-point methods):严格可行初始点通常难找;只要求 \(x^0,s^0>0\)。定义残差
路径跟踪方法(Path-Following Methods):显式限制迭代点在 \(\mathcal C\) 的邻域内并沿 \(\mathcal C\) 走向解,防止过于接近非负象限边界,保证每步方向至少有最低限度的进展。对偶度量 \(\mu\) 充当优劣度量,迫使 \(\mu_k\to0\)。两种邻域:
Algorithm 14.2(长步路径跟踪,Long-Step Path-Following):给定 \(\gamma\in(0,1)\),\(0<\sigma_{\min}<\sigma_{\max}<1\),\((x^0,\lambda^0,s^0)\in\mathcal N_{-\infty}(\gamma)\);对 \(k=0,1,\dots\):选 \(\sigma_k\in[\sigma_{\min},\sigma_{\max}]\);解 (14.12);取 \(\alpha_k\) 为 \([0,1]\) 中使 \((x^k(\alpha),\lambda^k(\alpha),s^k(\alpha))\in\mathcal N_{-\infty}(\gamma)\) 的最大值((14.19));更新。 图 14.2(\(n=2\)):横纵轴为 \(x_1s_1\)、\(x_2s_2\),中心路径为过原点 45° 线;在此几何下搜索方向变成曲线。\(\sigma_{\min}\) 保证每个方向起初离开 \(\mathcal N_{-\infty}(\gamma)\) 边界、进入其相对内部(小步改善中心性);较大 \(\alpha\) 又会走出邻域(线性化 (14.11) 对非线性系统 (14.9) 的误差随 \(\alpha\) 增大),但保证能走一个最小步长。14.4 节给出完整分析,说明原始–对偶理论无需深奥数学。该算法实践中相当高效,再加几项改动即成为真正有竞争力的方法。不可行变体:把 \(\mathcal N_{-\infty}(\gamma)\) 推广为允许违反可行性,要求 \(\|r_b\|,\|r_c\|\) 不超过 \(\mu\) 的常数倍,压低 \(\mu\) 即同时迫使残差趋零。
14.2 实用的原始–对偶算法(A Practical Primal–Dual Algorithm)(PDF p.421–426)
多数通用 LP 内点代码基于 Mehrotra 预测–校正算法(predictor–corrector),两大特点:(a) 在 Framework 14.1 的方向上加校正步(corrector),更紧密地跟踪通往解集的轨迹;(b) 自适应选择中心化参数 \(\sigma\)。 动机:把中心路径"平移"为从当前点 \((x,\lambda,s)\) 出发、终于解集 \(\Omega\) 的轨迹 \(\mathcal H=\{(\hat x(\tau),\hat\lambda(\tau),\hat s(\tau)):\tau\in[0,1)\}\),\(\tau=0\) 时为当前点,\(\tau\uparrow1\) 时极限属于 \(\Omega\)(图 14.3)。Framework 14.1 的算法可视为一阶方法:求轨迹的切线(预测步,predictor)并沿之线搜索。Mehrotra 进一步算出 \(\mathcal H\) 在当前点的曲率,得到二阶近似,曲率定义校正步;重用预测步的矩阵分解,边际成本低。自适应 \(\sigma\):先算仿射尺度方向(预测步)并评估其效果;若它能大幅降低 \(\mu\) 而不违反正性,说明不需多少中心化,取 \(\sigma\) 接近 0;否则取 \(\sigma\) 接近 1。最终方向由三部分组成:预测步(确定 \(\sigma_k\))、利用 \(\mathcal H\) 二阶信息的校正步、以 \(\sigma_k\) 代入 (14.15) 的中心化步;中心化与校正可合并计算,自适应中心化不增加每步代价。
具体计算:
- 预测步:在 (14.15) 中取 \(\sigma=0\):
\[\begin{bmatrix}0&A^T&I\\A&0&0\\S&0&X\end{bmatrix}\begin{bmatrix}\Delta x^{\rm aff}\\\Delta\lambda^{\rm aff}\\\Delta s^{\rm aff}\end{bmatrix}=\begin{bmatrix}-r_c\\-r_b\\-XSe\end{bmatrix}.\quad(14.20)\]
- 沿此方向不违反非负性的最大步长(上限 1):
\[\alpha^{\rm pri}_{\rm aff}=\min\Big(1,\min_{i:\Delta x_i^{\rm aff}<0}-\frac{x_i}{\Delta x_i^{\rm aff}}\Big),\quad\alpha^{\rm dual}_{\rm aff}=\min\Big(1,\min_{i:\Delta s_i^{\rm aff}<0}-\frac{s_i}{\Delta s_i^{\rm aff}}\Big),\quad(14.21)\]\[\mu_{\rm aff}=(x+\alpha^{\rm pri}_{\rm aff}\Delta x^{\rm aff})^T(s+\alpha^{\rm dual}_{\rm aff}\Delta s^{\rm aff})/n,\quad(14.22)\](原书 (14.22) 第二个因子也写作 \(\alpha^{\rm pri}_{\rm aff}\),按 Mehrotra 原意应为对偶步长。)
- 中心化参数 \(\sigma=(\mu_{\rm aff}/\mu)^3\):预测方向进展好时 \(\mu_{\rm aff}\ll\mu\),\(\sigma\) 小;反之 \(\sigma\) 大。
- 校正步右端为 \((0,0,-\Delta X^{\rm aff}\Delta S^{\rm aff}e)\)(补偿线性化忽略的二阶项 \(\Delta x_i\Delta s_i\)),中心化步右端为 \((0,0,\sigma\mu e)\)。三者右端相加,一次求解:
\[\begin{bmatrix}0&A^T&I\\A&0&0\\S&0&X\end{bmatrix}\begin{bmatrix}\Delta x\\\Delta\lambda\\\Delta s\end{bmatrix}=\begin{bmatrix}-r_c\\-r_b\\-XSe-\Delta X^{\rm aff}\Delta S^{\rm aff}e+\sigma\mu e\end{bmatrix}.\quad(14.23)\]
- 最大步长 \(\alpha^{\rm pri}_{\max}=\min(1,\min_{i:\Delta x_i<0}-x_i^k/\Delta x_i)\),\(\alpha^{\rm dual}_{\max}\) 同理((14.24));实际步长
\[\alpha^{\rm pri}_k=\min(1,\eta\alpha^{\rm pri}_{\max}),\quad\alpha^{\rm dual}_k=\min(1,\eta\alpha^{\rm dual}_{\max}),\quad(14.25)\]\(\eta\in[0.9,1.0)\),接近解时 \(\eta\to1\) 以加速渐近收敛。原始与对偶分别取步长。
Algorithm 14.3(Mehrotra 预测–校正):给定 \((x^0,\lambda^0,s^0)\),\((x^0,s^0)>0\);每步:解 (14.20) 得仿射方向;按 (14.21)(14.22) 算 \(\alpha^{\rm pri}_{\rm aff},\alpha^{\rm dual}_{\rm aff},\mu_{\rm aff}\);\(\sigma=(\mu_{\rm aff}/\mu)^3\);解 (14.23);按 (14.25) 算步长;\(x^{k+1}=x^k+\alpha_k^{\rm pri}\Delta x\),\((\lambda^{k+1},s^{k+1})=(\lambda^k,s^k)+\alpha_k^{\rm dual}(\Delta\lambda,\Delta s)\)。 注意:上述形式的 Mehrotra 算法没有收敛理论,甚至存在发散的例子;可加简单保护措施纳入已有收敛框架,但多数程序不实现,因为其实际表现很好。\(\eta\) 和初始点的选择细节见 Mehrotra 原文。
求解线性系统(Solving the Linear Systems):主要计算量在于解 (14.15)、(14.20)、(14.23);系数矩阵大而稀疏,但可改写为更紧凑的对称形式。因 \(x,s>0\),\(X,S\) 非奇异。从 (14.15) 消去 \(\Delta s\):
14.3 其他原始–对偶算法与推广(Other Primal–Dual Algorithms and Extensions)(PDF p.426–428)
其他路径跟踪方法:
- 短步路径跟踪(short-step):取保守的 \(\sigma\)(略小于 1),使单位步 \(\alpha=1\) 不离开限制性 \(\mathcal N_2\) 邻域 (14.16),进展慢。
- Mizuno–Todd–Ye 预测–校正法(与 Algorithm 14.3 不同):用两个嵌套的 \(\mathcal N_2\) 邻域。每隔一步为预测步:从内邻域出发沿仿射尺度方向(\(\sigma=0\))走到外邻域边界,两边界间隙足以大幅降低 \(\mu\);交替的校正步(\(\sigma=1\),\(\alpha=1\))把迭代拉回内邻域。\(\mu_k\) 超线性收敛到零(多数方法为线性)。
势函数下降方法(Potential-Reduction Methods):步的形式与路径跟踪相同,但不显式跟随中心路径;用对数势函数衡量 \(\mathcal F^o\) 中点的优劣,每步要求势函数固定量下降。原始–对偶势函数 \(\Phi\) 通常有两性质:
推广(Extensions):
- 单调线性互补问题(LCP):求 \(x,s\in\mathbb R^n\) 使
\[s=Mx+q,\quad(x,s)\ge0,\quad x^Ts=0,\quad(14.31)\]\(M\) 半正定。与 KKT (14.3) 形式相似。实例见 Cottle–Pang–Stone。
- 凸二次规划:
\[\min c^Tx+\tfrac12x^TGx\ \text{s.t.}\ Ax=b,\ x\ge0,\quad(14.32)\]\(G\) 对称半正定,KKT 条件类似 (14.3) 与 LCP;任何 LCP 可写成凸 QP,反之亦然(16.7 节讨论 QP 内点法)。这两类推广保留 LP 算法的收敛性与多项式复杂度。
- 非线性规划:把 KKT 写成类似 (14.3) 的形式(必要时加松弛把不等式化为简单界),对等式 KKT 条件用牛顿法并截短步长保持严格正性。凸问题的全局收敛不难证明,一般非线性情形见 17.2 节。
- 半定规划(semidefinite programming, SDP):部分变量为须半正定的对称矩阵,内点法很有效,应用于控制理论、组合优化等(Nesterov–Nemirovskii、Boyd 等、Vandenberghe–Boyd)。
14.4 Algorithm 14.2 的分析(Analysis of Algorithm 14.2)(PDF p.428–433)
从一个纯技术引理出发,几页内得到有力的定理:引理 14.1 → 引理 14.2(乘积向量的界)→ 定理 14.3(\(\alpha_k\) 下界与每步 \(\mu\) 的下降,蕴含全局收敛)→ 定理 14.4(\(O(n\log1/\epsilon)\) 次迭代达到 \(\mu_k<\epsilon\))。
引理 14.1:\(u,v\in\mathbb R^n\),\(u^Tv\ge0\),则 \(\|UVe\|\le2^{-3/2}\|u+v\|^2\),\(U=\mathrm{diag}(u)\),\(V=\mathrm{diag}(v)\)。 证明:对 \(\alpha\beta\ge0\) 有 \(\sqrt{|\alpha\beta|}\le\tfrac12|\alpha+\beta|\)((14.33),算术–几何平均)。划分 \(\mathcal P=\{i:u_iv_i\ge0\}\),\(\mathcal M=\{i:u_iv_i<0\}\),由 \(u^Tv\ge0\) 得 \(\sum_{\mathcal P}|u_iv_i|\ge\sum_{\mathcal M}|u_iv_i|\)((14.34))。于是 \(\|UVe\|=(\|[u_iv_i]_{\mathcal P}\|^2+\|[u_iv_i]_{\mathcal M}\|^2)^{1/2}\le(\|[\cdot]_{\mathcal P}\|_1^2+\|[\cdot]_{\mathcal M}\|_1^2)^{1/2}\le(2\|[u_iv_i]_{\mathcal P}\|_1^2)^{1/2}\le\sqrt2\|[\tfrac14(u_i+v_i)^2]_{\mathcal P}\|_1=2^{-3/2}\sum_{\mathcal P}(u_i+v_i)^2\le2^{-3/2}\|u+v\|^2\)。
引理 14.2:若 \((x,\lambda,s)\in\mathcal N_{-\infty}(\gamma)\),则 (14.12) 的步满足 \(\|\Delta X\Delta Se\|\le2^{-3/2}(1+1/\gamma)n\mu\)。(抽取文本中 \(\Delta\) 符号丢失,此处及证明中的乘积均指步向量的 \(\Delta X\Delta Se\)。) 证明:由 (14.12) 前两行得 \(\Delta x^T\Delta s=0\)((14.35),\(\Delta x\in\mathrm{Null}(A)\),\(\Delta s=-A^T\Delta\lambda\))。第三行乘 \((XS)^{-1/2}\),用 \(D=X^{1/2}S^{-1/2}\):
定理 14.3:给定 Algorithm 14.2 的 \(\gamma,\sigma_{\min},\sigma_{\max}\),存在与 \(n\) 无关的常数 \(\delta\) 使
定理 14.4(复杂度):给定 \(\epsilon>0\)、\(\gamma\in(0,1)\),初始点属于 \(\mathcal N_{-\infty}(\gamma)\) 且
第 14 章 注释与参考(PDF p.433–434)
详见 Wright 的专著《Primal-Dual Interior-Point Methods》。Karmarkar 方法源于寻找最坏情况优于单纯形法的 LP 算法;首个多项式算法(Khachiyan 椭球法)计算上令人失望,而 Karmarkar 方法提出时执行时间与当时单纯形代码相差不大(大问题上尤其如此)。Karmarkar 算法是纯原始算法:每步对原始可行集做投影变换,把当前点映到集合中心,再沿变换空间中的可行最速下降方向走一步,用对数势函数度量进展(见 Karmarkar 原文、Fletcher);其实际表现似乎不及最有效的原始–对偶方法。本章的路径跟踪与势函数方法也都具有多项式复杂度。 历史:许多 1984 年后研究的思想源于此前三部工作——Fiacco–McCormick 关于对数障碍函数的书(证明了中心路径存在等);McLinden 在非线性互补问题背景下对中心路径的进一步分析;Dikin 提出的**原始仿射尺度(primal affine-scaling)**内点法。原始–对偶方法的研究热潮始于 Megiddo 的开创性论文。Todd 对势函数下降法有出色综述(联系纯原始势函数法,包括 Karmarkar 原始算法)。复杂性理论入门见 Vavasis。LP 内点软件已广泛可得,多数基于 Algorithm 14.3 加 Gondzio 提出的"高阶校正"步;因比单纯形代码简单,部分可免费用于研究甚至商业用途。
第 14 章习题(PDF p.434–436)
- 14.1 LP \(\min x_1\) s.t. \(x_1+x_2=1\),\(x\ge0\):原始–对偶解 \(x^*=(0,1)\),\(\lambda^*=0\),\(s^*=(1,0)\);验证 \(F=0\) 有伪解 \(x=(1,0)\),\(\lambda=1\),\(s=(0,-1)\),与 LP 解无关(说明非负界的必要性)。
- 14.2 (i) \(\mathcal N_2(\theta_1)\subset\mathcal N_2(\theta_2)\)(\(\theta_1<\theta_2\)),\(\mathcal N_{-\infty}(\gamma_1)\subset\mathcal N_{-\infty}(\gamma_2)\)(\(\gamma_2\le\gamma_1\));(ii) \(\gamma\le1-\theta\) 时 \(\mathcal N_2(\theta)\subset\mathcal N_{-\infty}(\gamma)\)。
- 14.3 给定 \(\mathcal F^o\) 中一点,求使其属于 \(\mathcal N_{-\infty}(\gamma)\) 的 \(\gamma\) 范围。14.4 \(n=2\) 时构造不属于任何 \(\mathcal N_2(\theta)\) 的点。14.5 \(\mathcal N_{-\infty}(1)\) 与 \(\mathcal N_2(0)\) 都等于中心路径。
- 14.6 证 \(\Phi_\rho\) 满足 (14.29a)。14.7 (14.11) 系数矩阵非奇异 ⇔ \(A\) 行满秩。14.8 证 (14.35)。
- 14.9 \(AD^2A^T\) 对称正定 ⇔ \(A\) 行满秩;若 \(D\) 恰有 \(m\) 个正对角元其余为零是否仍成立。
- 14.10 对轨迹 \(\mathcal H\)(\(F(\hat x(\tau),\hat\lambda(\tau),\hat s(\tau))=(1-\tau)(r_c,r_b,XSe)\))求 \(\tau=0\) 处一、二、三阶导数方程,写出 Taylor 近似(预测–校正的理论基础)。
- 14.11 含自由变量 \(y\) 的 LP \(\min c^Tx+d^Ty\) s.t. \(A_1x+A_2y=b\),\(x\ge0\):写最优性条件、步方程、增广系统形式,解释为何不能化为对称正定的正规方程形式。
- 14.12 用 Matlab 实现 Algorithm 14.3(\(\eta=0.99\)),用随机 \(A\) 构造已知解的 LP(\(x\) 前 \(m\) 分量随机正、其余 0;\(s\) 相反;\(\lambda\) 随机;\(c=A^T\lambda+s\),\(b=Ax\)),初始 \(x^0,s^0\) 取大正数。
第 14 章本章要点
- 原始–对偶内点法对 KKT 等式 \(F(x,\lambda,s)=0\) 用牛顿法,并保持 \((x,s)>0\) 严格成立以避开伪解。
- 中心路径:\(x_is_i=\tau\) 的严格可行点曲线;对偶度量 \(\mu=x^Ts/n\);中心化参数 \(\sigma\) 在仿射尺度(\(\sigma=0\))与中心化(\(\sigma=1\))之间权衡。
- 路径跟踪:迭代保持在 \(\mathcal N_2(\theta)\) 或 \(\mathcal N_{-\infty}(\gamma)\) 邻域内;长步法每步 \(\mu\) 至少下降因子 \(1-\delta/n\),\(O(n\log1/\epsilon)\) 次迭代——多项式复杂度。
- 实用算法:Mehrotra 预测–校正(仿射预测 → \(\sigma=(\mu_{\rm aff}/\mu)^3\) → 二阶校正 + 中心化合并一次求解 → 原始/对偶分别取步长 \(\eta\alpha_{\max}\)),无收敛理论但实践中极好。
- 主要计算是解正规方程 \(AD^2A^T\Delta\lambda=\cdots\)(稀疏 Cholesky),后期病态;增广系统处理稠密列与自由变量更好。
- 推广到 LCP、凸 QP、非线性规划与半定规划。
第 14 章与量化交易的关联
- 组合优化求解器的内核:MOSEK、Gurobi、CPLEX 的 barrier 算法,以及 CVXPY 背后的 ECOS、Clarabel、SCS(后者为一阶法)都是本章方法的推广;均值-方差(凸 QP)、带跟踪误差/风险预算的二阶锥规划(SOCP)、协方差矩阵最近相关阵修复(SDP)都是内点法的典型应用。理解中心路径、\(\mu\)、原始/对偶残差有助于读懂求解器日志、设置容差和诊断"数值问题/接近不可行"警告。
- 大规模问题的规模化:数千资产、多期、带大量线性约束的组合优化中,内点法迭代次数几乎与规模无关(通常 20–60 次),而单纯形法可能随规模迅速增长;但内点法给出的是"内部"近似解,需交叉(crossover)才能得到顶点/基解——这对需要精确稀疏持仓的场景有影响。
- 正规方程 \(AD^2A^T\) 与稠密列:组合问题中预算约束 \(\mathbf 1^Tw=1\) 或因子暴露约束常是稠密行/列,会导致正规方程填充;了解增广系统形式有助于在自写求解器(如定制 ADMM/内点)时规避性能陷阱。
- LCP:美式期权离散化后的定价问题(障碍问题)可写成线性互补问题,可用本章内点思想或投影 SOR 求解。
- 对偶间隙作停止准则:内点法天然给出对偶界,可量化"当前组合离最优还有多远",适合需要时效的盘中再优化。
第 14 章推荐习题
- 14.1(伪解,理解为何必须保持内点);
- 14.5、14.2(中心路径与邻域的关系);
- 14.9、14.11(正规方程与增广系统的适用条件,自由变量处理);
- 14.12(亲手实现 Mehrotra 算法,最有实践价值);
- 14.10(预测–校正的 Taylor 展开推导)。
第 15 章 非线性约束优化算法基础(Fundamentals of Algorithms for Nonlinear Constrained Optimization)(PDF p.437–444,续见下一块)
章引言(PDF p.438–441):开始讨论一般约束优化
- 线性规划:\(f\)、\(c_i\) 全线性(第 13、14 章);
- 二次规划:约束线性、目标二次(第 16 章);
- 非线性规划:至少部分约束是一般非线性函数;
- 线性约束优化:所有约束线性;
- 界约束优化:约束只有 \(x_i\ge l_i\) 或 \(x_i\le u_i\);
- 凸规划:\(f\) 凸,等式约束线性,不等式约束 \(c_i\) 凹。 这些类别既不互斥也不穷尽,且可细分(如凸二次规划是目标凸的 QP 子类);更细的分类对算法有意义,例如凸问题更容易选价值函数。 约束优化算法都是迭代的:产生趋于解 \(x^*\) 的猜测序列,也可能产生 Lagrange 乘子的估计序列;利用目标、约束及其导数信息(可能结合历史信息)决定下一步;找到近似解或无法继续进展时终止。本书只研究求局部极小点的算法,全局极小不在范围内。
问题的初步研究(Initial Study of a Problem):求解前先研究问题能否简化。有时无需计算机即可求解,如检查约束发现可行域为空或目标在可行域上无界(习题 15.1)。也可能猜测哪些不等式在解处活跃,把 KKT 条件化为可直接求解的方程组(如第 12 章若干例子);但这很少实用——即使能识别活跃约束(通常是约束优化算法面临的最难问题),仍需数值求解方程组,而第 11 章表明非线性方程算法不能保证从任意起点找到解。所以需要直接处理 (15.1) 的非线性优化算法。 确定问题类别;若含离散变量(如 0/1 二值变量),本书技术不适用,需用离散优化算法(Wolsey;Nemhauser–Wolsey)。 硬约束与软约束:从算法角度,**硬约束(hard constraints)**是为使 (15.1) 中函数有意义而必须满足的约束,某些函数在不可行点甚至无定义(例:变量须为正,因目标中需要其平方根;或所有变量之和须为零以满足守恒律)。软约束(soft constraints)有时被建模者改写为无约束问题——在目标中加入含约束的罚项;下一章将看到罚方法通常引入病态,是否有害取决于所用无约束算法。用户需权衡显式处理约束还是用罚方法。 必须在所有迭代点满足硬约束时须用可行算法(feasible algorithms):选一个满足硬约束的初始点,并保持新迭代点对这些约束可行。可行算法通常比允许不可行迭代的算法更慢、更贵(不能抄穿越不可行区域的近路),但优点是可直接用 \(f\) 判断每个点的优劣,无需引入考虑约束违反的复杂价值函数。
15.1 优化算法的分类(Categorizing Optimization Algorithms)(PDF p.441–443)
非线性优化算法没有"标准分类",本书其余章节分组如下: I. 第 17 章:罚、障碍、增广 Lagrange 方法,以及序列线性约束方法。
- 二次罚函数:只有等式约束时
\[f(x)+\frac1{2\mu}\sum_{i\in\mathcal E}c_i^2(x),\quad(15.2)\]\(\mu>0\) 为罚参数;对一系列递减的 \(\mu\)(即递增的罚权重 \(1/\mu\);原文写"increasing values of \(\mu\)",与 \(1/(2\mu)\) 的写法不一致,应理解为罚权重增加)极小化该无约束函数,直到足够精确地识别约束问题的解。
- 精确罚函数:可能一次无约束极小化即可解 (15.1)。等式约束情形 \(f(x)+\frac1\mu\sum_{i\in\mathcal E}|c_i(x)|\),\(\mu\) 足够小(但为正)。精确罚函数通常不可微,极小化它需解一系列子问题。
- 障碍方法:在目标中加入当 \(x\) 处于可行集内部时不显著、但 \(x\) 接近边界时趋于无穷的项。只有不等式约束时,对数障碍函数 \(f(x)-\mu\sum_{i\in\mathcal I}\log c_i(x)\),\(\mu>0\) 为障碍参数;在一定条件下其极小点随 \(\mu\downarrow0\) 趋于原问题的解;对递减的 \(\mu\) 序列求近似极小点。(原文称障碍项"approach zero",实际上 \(-\log c_i\to+\infty\)。)
- 增广 Lagrange 方法:结合 Lagrange 函数 (12.28) 与二次罚函数 (15.2)。只有等式约束时
\[\mathcal L_A(x,\lambda;\mu)=f(x)-\sum_{i\in\mathcal E}\lambda_ic_i(x)+\frac1{2\mu}\sum_{i\in\mathcal E}c_i^2(x).\]固定 \(\lambda\)(最优乘子估计)与 \(\mu>0\),近似极小化 \(\mathcal L_A\) 得新 \(x\),用之更新 \(\lambda\),可能减小 \(\mu\),重复。可避免罚函数和障碍函数极小化的某些数值困难(病态)。
- 序列线性约束方法(sequential linearly constrained):每步在约束线性化下极小化某个 Lagrange 函数,主要用于大规模问题。
II. 第 18 章:序列二次规划(SQP)。 每个迭代点用二次子问题建模 (15.1),搜索方向为子问题的解。等式约束情形,在 \((x_k,\lambda_k)\) 处
\[\min_p\ \tfrac12p^TW_kp+\nabla f_k^Tp\quad(15.3a)\qquad\text{s.t. }A_kp+c_k=0,\quad(15.3b)\]目标近似 Lagrange 函数(\(W_k\) 为其 Hessian 或近似),约束为线性化约束;沿该方向搜索直到某价值函数下降。SQP 实践中很有效,是求解大小规模约束优化问题的一些最佳软件的基础;通常比其他方法需要更少的函数求值,代价是每步要解相对复杂的二次子问题。 III. 第 16 章:二次规划算法。 因其重要且算法可针对其特点定制,单独成章:有效集方法与内点法;有效集 QP 方法是上述 SQP 的基础。 II、III 类算法使用约束消元技术,下面先讨论这一背景;本章随后讨论价值函数(SQP 等算法的重要组件)。 学习提示:以下概念属背景材料,读者可只浏览后两节,在学习第 17、18 章时按需回看。
15.2 变量消元(Elimination of Variables)(PDF p.443–444,续见下一块)
自然思路:消去约束得到无约束问题,或至少消去部分约束得到更简单的问题。但消元须谨慎,可能改变问题或引入病态。 安全的例子:\(\min f(x_1,x_2,x_3,x_4)\) s.t. \(x_1+x_3^2-x_4x_3=0\),\(-x_2+x_4+x_3^2=0\)(原文第一个约束印为 \(x_4x_5\),按代入式应为 \(x_4x_3\)),可令 \(x_1=x_4x_3-x_3^2\),\(x_2=x_4+x_3^2\),极小化二元无约束函数 \(h(x_3,x_4)=f(x_4x_3-x_3^2,\ x_4+x_3^2,\ x_3,\ x_4)\),用前面章节的任何无约束算法。 例 15.1(非线性消元的危险,Fletcher):\(\min x^2+y^2\) s.t. \((x-1)^3=y^2\),解为 \((1,0)\)(图 15.1)。消去 \(y\) 得 \(h(x)=x^2+(x-1)^3\),\(x\to-\infty\) 时 \(h\to-\infty\),盲目变换会误以为问题无界;但约束 \((x-1)^3=y^2\ge0\) 隐含界 \(x\ge1\),且在解处活跃,消元时必须显式加入此界(图 15.1 的图注即强调:消元后应把错误纳入的那一支从可行集中排除)。 因此非线性消元可能产生难以追踪的错误,多数优化算法不用;而是先线性化约束,再对简化问题做消元。接下来(下一块)系统地介绍线性约束的消元过程。(续见下一块)