量化交易中文教材

第 17 章 罚函数、障碍函数与增广拉格朗日方法

学习目标

读完本章,你应当能够:

  1. 写出二次罚函数 \(Q(x;\mu)=f+\frac1{2\mu}\|c\|^2\),证明其极小点在 \(\mu\to0\) 时收敛到 KKT 点,并用 \(-c_i(x_k)/\mu_k\) 估计乘子;说明它的 Hessian 为什么会严重病态。
  2. 写出对数障碍函数 \(P(x;\mu)=f-\mu\sum\log c_i\),理解中心路径、乘子估计 \(\mu/c_i\) 与"扰动互补" \(\lambda_ic_i=\mu\),并把它与第 14、16 章的原始–对偶内点法联系起来。
  3. 说明精确罚函数(\(\ell_1\))的优缺点。
  4. 推导增广拉格朗日函数与乘子更新 \(\lambda\leftarrow\lambda-c/\mu\),理解为什么它不需要 \(\mu\to0\)(定理 17.5、17.6),会用 \(\psi\) 函数或松弛 + 界约束处理不等式。
  5. 了解 LANCELOT 的实用乘子法和序列线性约束方法(MINOS)。
  6. 在量化场景中正确使用罚项(软约束)、障碍项(风险平价)和增广拉格朗日 / ADMM。

读前导读

这一章在解决什么问题。 有约束的优化很难直接解,而无约束优化(第 2–9 章的牛顿法、BFGS)已经很成熟。本章的思路是"把约束搬进目标函数":违反约束就罚钱,或者靠近边界就收高额"过路费",然后反复解无约束问题。你在组合管理里其实天天这么做:把跟踪误差、换手成本写成 \(\lambda\cdot\text{TE}^2\) 加进均值–方差目标,这就是二次罚函数。本章回答三个实际问题:罚系数为什么不能无脑调大(会病态);为什么"罚 + 乘子校正"(增广拉格朗日)能用不大的罚系数就精确满足约束;对数障碍为什么能让迭代点始终可行,以及它和 MOSEK、Gurobi 的"barrier"算法是同一个东西。

拉格朗日乘子在这里始终扮演影子价格的角色:约束放松一单位,最优目标改善多少。罚函数法和障碍法都会"顺带"给出这个影子价格的估计(\(-c_i/\mu\) 和 \(\mu/c_i\)),增广拉格朗日法则是主动去修正这个价格。理解了这一点,全章就是一条线。

需要先想起来的数学。

  • 梯度与 Hessian、链式法则。 \(\nabla f\) 是各偏导组成的向量,\(\nabla^2 f\) 是二阶偏导组成的矩阵。本章大量用链式法则:\(\nabla(c(x)^2)=2c(x)\nabla c(x)\),\(\nabla\log c(x)=\nabla c(x)/c(x)\)。例如 \(c(x)=x_1^2+x_2^2-2\),则 \(\nabla c=(2x_1,2x_2)\),在 \((-1,-1)\) 处为 \((-2,-2)\)。见 第 00 册第 05 章 多元微积分与优化。
  • 拉格朗日乘子与 KKT 条件。 本书沿用原书记号 \(\mathcal L(x,\lambda)=f(x)-\sum_i\lambda_ic_i(x)\)(注意是减号),驻点条件为 \(\nabla f=\sum_i\lambda_i\nabla c_i\)。不等式 \(c_i\ge0\) 要求 \(\lambda_i\ge0\) 且 \(\lambda_ic_i=0\)(互补松弛:约束不绑定则影子价格为零)。\(\mathcal E\)、\(\mathcal I\) 分别是等式和不等式约束的下标集合。LICQ 指有效约束的梯度线性无关。详见第 12 章和 第 00 册第 05 章。
  • 矩阵的特征值与条件数。 对称正定矩阵的条件数 = 最大特征值 / 最小特征值,衡量"各方向的曲率差多少倍"。例如 \(\mathrm{diag}(1,1000)\) 的条件数是 1000,等高线是极扁的椭圆,梯度类方法会来回振荡。\(A(x)\) 是约束的雅可比矩阵(第 \(i\) 行是 \(\nabla c_i(x)^T\)),\(A^TA\) 半正定、秩不超过约束个数。\(\mathrm{Null}\,A=\{w:Aw=0\}\) 是零空间,即"沿着它移动、线性化约束不变"的方向。见 第 00 册第 06 章 线性代数速成。
  • 大 O 记号与极限记号。 \(O(1/\mu)\) 表示"量级与 \(1/\mu\) 相同";\(\mu_k\downarrow0\) 表示 \(\mu_k\) 单调递减趋于 0;"极限点"指序列中某个子序列收敛到的点。见 第 00 册第 07 章 和 第 01 章。

怎么读这一章。 核心必读是 17.2(二次罚,尤其 17.2.2 的乘子估计和 17.2.3 的病态)、17.3.1–17.3.4(对数障碍与扰动互补)和 17.5.1(增广拉格朗日的动机与乘子更新)。这三处把握住,再看 17.7 的两段代码输出,就能把理论和数字对上。第一次可以只看结论的部分:17.3.5 的原始–对偶矩阵 (17.43)(知道它是 LP/QP 内点法的通式即可)、17.5.2 的证明细节 (17.62)、17.5.4 LANCELOT 和 17.6 MINOS。建议顺序:17.1 → 17.2 → 17.5.1 →(回头)17.3 → 17.7 → 其余。先读 17.5.1 是因为它直接修补 17.2 的缺陷,两节连读最顺。


17.1 核心思路

本章的方法有一个共同点:用一列无约束(或只有简单界约束的)子问题代替原约束问题。三种主要方法:

  • 二次罚函数法(quadratic penalty method):把约束违反量的平方乘一个系数加到目标上。简单直观、广泛使用,但有病态问题。
  • 对数障碍法(log-barrier method):用对数项阻止可行迭代点靠近边界。它是 LP、凸 QP、半定规划的原始与原始–对偶内点法的基础。
  • 增广拉格朗日法 / 乘子法(augmented Lagrangian / method of multipliers):在罚函数中显式加入乘子估计,避免了二次罚函数固有的病态。

另有精确罚函数(一次极小化代替一列问题,但难以极小化)和序列线性约束方法(大规模问题的实用方法)。

量化映射:把跟踪误差、换手、因子暴露偏离写成罚项加进目标,是二次罚;风险平价 \(\min\frac12w^T\Sigma w-\sum b_i\log w_i\) 是障碍函数;ADMM(OSQP、大规模 Lasso 因子选择、分布式组合优化)是增广拉格朗日的分块版本;商业求解器(MOSEK、IPOPT、Gurobi barrier)的内核是原始–对偶障碍法。


17.2 二次罚函数法

17.2.1 定义

罚函数 = 原目标 + 每个约束一个罚项(违反时为正、满足时为零),罚项乘一个正系数;系数越来越大,极小点越来越接近可行域。由于罚函数的极小点通常不可行、只在极限下趋于可行,这类方法又叫外点罚方法(exterior penalty methods)。

等式约束问题

\[ \min_x f(x)\quad\text{s.t.}\quad c_i(x)=0,\ i\in\mathcal E\tag{17.1} \]

的二次罚函数为

\[ Q(x;\mu)=f(x)+\frac1{2\mu}\sum_{i\in\mathcal E}c_i^2(x),\tag{17.2} \]

\(\mu>0\) 为罚参数,令 \(\mu_k\downarrow0\)。罚项是光滑的(\(c_i^2\) 与 \(c_i\) 同阶可微),可以用任何无约束方法;上一个 \(\mu\) 的近似极小点给下一个 \(\mu\) 提供好的起点。

白话解释:\(\frac1{2\mu}\) 就是"罚款单价",\(\mu\) 越小单价越高。约束违反 \(c_i\) 取平方是为了两个方向都罚、而且在 \(c_i=0\) 处光滑。为什么不一次把 \(\mu\) 设得极小?因为单价一高,目标函数在约束附近形成又窄又陡的"峡谷",优化器很难走(17.2.3 的病态)。所以做法是从温和的罚款开始,逐步加码,每次用上一轮的解热启动——和你做敏感性分析时逐步收紧约束、每次从上一组权重出发是同一个习惯。

例 17.1:\(\min x_1+x_2\) s.t. \(x_1^2+x_2^2-2=0\),解 \(x^*=(-1,-1)\),\(\lambda^*=-\frac12\)。

\[ Q(x;\mu)=x_1+x_2+\frac1{2\mu}(x_1^2+x_2^2-2)^2.\tag{17.4} \]

\(\mu=1\) 时极小点约在 \((-1.1,-1.1)\),另在 \((0.3,0.3)\) 附近有局部极大;\(\mu=0.1\) 时不在圆上的点受重罚,出现一道明显的低值"槽",极小点更接近 \((-1,-1)\)。

含不等式:

\[ Q(x;\mu)=f(x)+\frac1{2\mu}\sum_{i\in\mathcal E}c_i^2(x)+\frac1{2\mu}\sum_{i\in\mathcal I}([c_i(x)]^-)^2,\qquad[y]^-=\max(-y,0).\tag{17.5} \]

此时 \(Q\) 一般只有一阶连续导数:例如 \(x_1\ge0\) 对应 \(\min(0,x_1)^2\),在 \(x_1=0\) 处二阶导不连续(习题 17.1),这会影响牛顿类方法。

框架 17.1(二次罚函数)

给定 μ0 > 0,容差 τ0 > 0,起点 x0^s;
for k = 0,1,2,...
    从 xk^s 出发求 Q(·;μk) 的近似极小点 xk,当 ‖∇Q(x;μk)‖ ≤ τk 时停止;
    若满足最终收敛测试:STOP,返回 xk;
    选新罚参数 μ_{k+1} ∈ (0, μk);选新起点 x_{k+1}^s;
end

\(\mu\) 的减小速度可以自适应:子问题难时温和减小(如 \(\mu_{k+1}=0.7\mu_k\)),容易时激进减小(如 \(0.1\mu_k\))。收敛理论只要求 \(\tau_k\to0\)。

17.2.2 收敛性

定理 17.1:若每个 \(x_k\) 是 \(Q(x;\mu_k)\) 的精确全局极小点且 \(\mu_k\downarrow0\),则 \(\{x_k\}\) 的每个极限点都是 (17.1) 的全局解。

证明:设 \(\bar x\) 为全局解(\(c(\bar x)=0\))。由 \(Q(x_k;\mu_k)\le Q(\bar x;\mu_k)=f(\bar x)\),

\[ \sum c_i^2(x_k)\le2\mu_k[f(\bar x)-f(x_k)].\tag{17.7} \]

沿收敛子列取极限,右端趋于零,故极限点 \(x^*\) 可行;又 \(f(x_k)\le f(\bar x)\),取极限 \(f(x^*)\le f(\bar x)\),所以 \(x^*\) 是全局解。∎(实际中每个子问题求全局极小很难,所以更有用的是下面的结果。)

定理 17.2:若 \(\tau_k\to0\)、\(\mu_k\downarrow0\),则对 \(\{x_k\}\) 的任一使 \(\nabla c_i(x^*)\) 线性无关的极限点 \(x^*\),\(x^*\) 是 (17.1) 的 KKT 点,且沿相应子列

\[ \lim_{k\in\mathcal K}\ -\frac{c_i(x_k)}{\mu_k}=\lambda_i^*,\quad i\in\mathcal E.\tag{17.8} \]

证明要点:

\[ \nabla_xQ(x_k;\mu_k)=\nabla f(x_k)+\sum_i\frac{c_i(x_k)}{\mu_k}\nabla c_i(x_k).\tag{17.9} \]

由 \(\|\nabla_xQ\|\le\tau_k\) 得 \(\|\sum c_i(x_k)\nabla c_i(x_k)\|\le\mu_k[\tau_k+\|\nabla f(x_k)\|]\),取极限并由线性无关得 \(c(x^*)=0\)。再记 \(\lambda^k=-c(x_k)/\mu_k\),(17.9) 变为 \(A(x_k)^T\lambda^k=\nabla f(x_k)-\nabla_xQ\),用最小二乘解出 \(\lambda^k\) 并取极限,得到 \(\nabla f(x^*)-A(x^*)^T\lambda^*=0\)。∎

推导拆解:为什么 \(-c_i/\mu\) 恰好是乘子?把 (17.9) 和 KKT 驻点条件并排写:

  • 罚函数极小点:\(\nabla f(x_k)+\sum_i\dfrac{c_i(x_k)}{\mu_k}\nabla c_i(x_k)\approx0\)(梯度近似为零,因为 \(\|\nabla_xQ\|\le\tau_k\))。
  • KKT 点:\(\nabla f(x^*)-\sum_i\lambda_i^*\nabla c_i(x^*)=0\)。

两式形状完全一样,只是第一式里 \(\nabla c_i\) 前面的系数是 \(c_i/\mu_k\),第二式是 \(-\lambda_i^*\)。对上号就得到 \(\lambda_i^*\approx-c_i(x_k)/\mu_k\)。(17.9) 本身只是对 \(\frac1{2\mu}c_i^2\) 用链式法则:\(\nabla\big(\tfrac1{2\mu}c_i^2\big)=\tfrac1{2\mu}\cdot2c_i\cdot\nabla c_i\)。"线性无关"的作用是保证 \(A^T\lambda=\) 某向量 有唯一解,从而 \(\lambda^k\) 的极限是确定的。

用例 17.1 验算:\(\lambda^*=-\tfrac12\),所以罚函数解满足 \(c(x_k)\approx\tfrac12\mu_k>0\),即解落在圆的外侧一点(17.7.2 输出中 \(\mu=1\) 时 \(c=0.45\))。直觉上,目标想让 \(x_1+x_2\) 更小,于是"花一点罚款"往外多走了一段。

两点重要推论:

  • \(-c_i(x_k)/\mu_k\) 是乘子估计。换句话说,罚函数极小点处的约束违反量满足 \(c_i(x_k)\approx-\mu_k\lambda_i^*\):违反量与 \(\mu_k\) 成正比,只有 \(\mu\to0\) 才消失。这是 17.5 节增广拉格朗日法的出发点。
  • 约束梯度线性相关时,极限点可能不可行;原问题不可行时,二次罚方法常收敛到 \(\|c(x)\|^2\) 的驻点。

17.2.3 病态

罚函数的 Hessian 为

\[ \nabla^2_{xx}Q(x;\mu_k)=\nabla^2f(x)+\sum_i\frac{c_i(x)}{\mu_k}\nabla^2c_i(x)+\frac1{\mu_k}A(x)^TA(x).\tag{17.15} \]

在极小点附近,由 (17.8),

\[ \nabla^2_{xx}Q(x;\mu_k)\approx\nabla^2_{xx}\mathcal L(x,\lambda^*)+\frac1{\mu_k}A(x)^TA(x),\tag{17.16} \]

即"与 \(\mu\) 无关的 Lagrange Hessian"加上"秩为 \(|\mathcal E|\)、非零特征值为 \(O(1/\mu_k)\) 的矩阵"。当 \(|\mathcal E|<n\) 时,一部分特征值是 \(O(1)\),另一部分是 \(O(1/\mu_k)\),条件数按 \(1/\mu_k\) 增长。后果:

白话解释:(17.15) 到 (17.16) 的一步是把 \(c_i/\mu_k\) 换成 \(-\lambda_i^*\)(用 (17.8)),于是前两项合成 \(\nabla^2f-\sum\lambda_i^*\nabla^2c_i=\nabla^2_{xx}\mathcal L\)。关键在第三项:\(A^TA/\mu\) 只在"离开约束面"的方向上有曲率,在"沿着约束面"的方向上为零。以例 17.1 为例,在 \((-1,-1)\) 附近,沿圆周切线方向曲率是 \(O(1)\),沿半径方向曲率约 \(8/\mu\)。\(\mu=10^{-4}\) 时两者差 \(10^4\)–\(10^5\) 倍,正是 17.7.2 输出的条件数。

这就像一个收益率曲线模型里,短端参数的敏感度比长端大上万倍:数值上一点点舍入误差就会被放大,梯度方法会在峡谷两壁之间来回弹,沿谷底前进得极慢。

  • 拟牛顿法、共轭梯度法在病态问题上表现很差;
  • 牛顿法对病态不太敏感,但 (1) 牛顿方程本身病态;(2) 二次 Taylor 模型只在极小点很小的邻域内准确——等高线呈"香蕉形"而非椭圆,牛顿步可能进展缓慢。

牛顿方程的稳定改写:引入辅助向量 \(\zeta=A(x)p/\mu_k\),

\[ \begin{bmatrix}\nabla^2f(x)+\sum_i\frac{c_i(x)}{\mu_k}\nabla^2c_i(x)&A(x)^T\\A(x)&-\mu_kI\end{bmatrix} \begin{bmatrix}p\\\zeta\end{bmatrix}=\begin{bmatrix}-\nabla_xQ(x;\mu)\\0\end{bmatrix},\tag{17.18} \]

它与原牛顿方程有相同的 \(p\),但 \(\mu_k\to0\) 时系数矩阵趋于一个条件良好的 KKT 型矩阵。

量化提示:软约束的罚系数。把跟踪误差、换手或因子暴露写成 \(\frac{\rho}2\|c(w)\|^2\) 加进目标时,\(\rho=1/\mu\) 越大越接近硬约束,但 (17.16) 说明 Hessian 条件数按 \(\rho\) 增长,求解越慢越不稳;而且由 (17.8),约束违反量约为 \(\lambda^*/\rho\)——永远不会精确为零。若需要精确满足,应改用增广拉格朗日(罚系数适中 + 乘子更新)或直接用约束求解器,而不是把罚系数调到极大。


17.3 对数障碍法

17.3.1 定义与例子

对只有不等式的问题

\[ \min_x f(x)\quad\text{s.t.}\quad c_i(x)\ge0,\ i\in\mathcal I,\tag{17.19} \]

记严格可行域 \(\mathcal F^o=\{x:c_i(x)>0,\ \forall i\}\),设非空。障碍函数在 \(\mathcal F^o\) 内光滑、在 \(\mathcal F^o\) 外为无穷、接近边界时趋于 \(+\infty\)。最重要的是对数障碍,组合函数为

\[ P(x;\mu)=f(x)-\mu\sum_{i\in\mathcal I}\log c_i(x),\tag{17.22} \]

\(\mu>0\) 为障碍参数,极小点记为 \(x(\mu)\)。与罚函数相反,障碍法的迭代点始终严格可行,所以属于内点法。

例 17.2:\(\min x\) s.t. \(x\ge0\),\(1-x\ge0\),\(P(x;\mu)=x-\mu\log x-\mu\log(1-x)\)。\(\mu\) 小时 \(P\) 在大部分可行集上接近 \(f\),只在两端窄窄的"边界层"中趋于无穷;\(x(\mu)\to x^*=0\)。

例 17.3:

\[ \min(x_1+0.5)^2+(x_2-0.5)^2\quad\text{s.t.}\quad x_1,x_2\in[0,1],\tag{17.23} \]
\[ P(x;\mu)=(x_1+0.5)^2+(x_2-0.5)^2-\mu[\log x_1+\log(1-x_1)+\log x_2+\log(1-x_2)].\tag{17.24} \]

解 \(x^*=(0,0.5)\),约束 \(x_1\ge0\) 有效。\(\mu=1,0.1\) 时极小点周围的等高线接近椭圆,多数无约束算法都能用;\(\mu=0.01\) 时极小点被推进边界层,等高线被拉长且不再是椭圆(靠边界一侧几乎是直线)。拉长意味着尺度差,拟牛顿法、最速下降法、CG 都表现差;非椭圆意味着二次模型不能很好刻画真实函数,牛顿法也只在 \(x(\mu)\) 的小邻域内快速收敛。

17.3.2 梯度、Hessian 与病态

\[ \nabla_xP=\nabla f(x)-\sum_i\frac{\mu}{c_i(x)}\nabla c_i(x),\tag{17.25} \]
\[ \nabla^2_{xx}P=\nabla^2f(x)-\sum_i\frac{\mu}{c_i(x)}\nabla^2c_i(x)+\sum_i\frac{\mu}{c_i^2(x)}\nabla c_i(x)\nabla c_i(x)^T.\tag{17.26} \]

推导拆解:两式都是对 \(-\mu\log c_i(x)\) 逐项求导。

  • 一阶:\(\log\) 的导数是 \(1/c\),再乘内层导数 \(\nabla c_i\)(链式法则),得 \(\nabla(-\mu\log c_i)=-\dfrac{\mu}{c_i}\nabla c_i\),即 (17.25)。
  • 二阶:对 \(-\dfrac{\mu}{c_i}\nabla c_i\) 再求导用乘积法则。对 \(\nabla c_i\) 求导得 \(-\dfrac{\mu}{c_i}\nabla^2c_i\);对 \(\dfrac1{c_i}\) 求导得 \(-\dfrac{1}{c_i^2}\nabla c_i\),与前面的负号相乘变正,得 \(+\dfrac{\mu}{c_i^2}\nabla c_i\nabla c_i^T\)。这就是 (17.26)。
  • 代入 \(\lambda_i^*\approx\mu/c_i\) 后,\(\dfrac{\mu}{c_i^2}=\dfrac{(\mu/c_i)^2}{\mu}\approx\dfrac{(\lambda_i^*)^2}{\mu}\),得到 (17.27)。

一维小例子:\(\min x\) s.t. \(x\ge0\),\(P=x-\mu\log x\),\(P'=1-\mu/x=0\) 得 \(x(\mu)=\mu\),\(P''=\mu/x^2=1/\mu\)。障碍越弱(\(\mu\) 越小),解越贴边界,曲率越陡。

在 \(x(\mu)\) 附近 \(\lambda_i^*\approx\mu/c_i(x)\),代入得

\[ \nabla^2_{xx}P(x;\mu)\approx\nabla^2_{xx}\mathcal L(x,\lambda^*)+\sum_i\frac{(\lambda_i^*)^2}{\mu}\nabla c_i(x)\nabla c_i(x)^T,\tag{17.27} \]

与二次罚的 (17.16) 结构相同:若解处有 \(t\) 个有效不等式(\(0<t<n\)),则 \(t\) 个特征值为 \(O(1/\mu)\),其余为 \(O(1)\),越来越病态。不过朴素的误差分析过于悲观:由于梯度和 Hessian 的特殊结构,直接消元算出的牛顿方向在 \(\mu\) 很小之前仍是好方向。

**框架 17.2(对数障碍)**与框架 17.1 相同,把 \(Q\) 换成 \(P\)。由于病态,牛顿法几乎是唯一有效的内层方法。为让牛顿法几步内进入收敛域:

  • 用保持 \(x\in\mathcal F^o\) 的线搜索(回溯时先保证 \(c_i(x+\alpha p)>0\))或信赖域;
  • 用好的起点:沿以往极小点外推,或沿中心路径的切向预测:
\[ \nabla^2_{xx}P(x;\mu)\dot x-\sum_{i\in\mathcal I}\frac1{c_i(x)}\nabla c_i(x)=0,\qquad x_k^s=x_{k-1}+(\mu_k-\mu_{k-1})\dot x;\tag{17.29–17.30} \]
  • 子问题不难且起点可靠时,激进地取 \(\mu_{k+1}=0.1\mu_k\) 或 \(0.2\mu_k\)。

对数障碍法由 Frisch 在 1950 年代提出、Fiacco 与 McCormick 在 60 年代系统研究(SUMT),因为与原始–对偶内点法的联系在 1984 年 Karmarkar 之后重新受到关注。

17.3.3 收敛性质

定理 17.3(凸规划):设 \(f\) 与 \(-c_i\) 均凸,\(\mathcal F^o\) 非空,解集 \(\mathcal M\) 非空有界,\(\mu_k\downarrow0\)。则 (i) 对任意 \(\mu>0\),\(P(x;\mu)\) 在 \(\mathcal F^o\) 上凸且取得极小点,任何局部极小都是全局极小;(ii) 极小点序列有收敛子列,极限点都在 \(\mathcal M\) 中;(iii) \(f(x(\mu_k))\to f^*\),\(P(x(\mu_k);\mu_k)\to f^*\)。

\(\mathcal M\) 为空(\(\min x\) s.t. \(-x\ge0\))或无界(\(\min x_1\) s.t. \(x_1\ge0\))时结论一般不成立。非凸问题中结论是局部的,且可能出现"障碍函数无下界"(习题 17.2)或"极小点序列收敛到非解点"(习题 17.3)。

17.3.4 与 KKT 条件的关系:扰动互补

在 \(x(\mu)\) 处 \(\nabla_xP=0\):

\[ \nabla f(x(\mu))-\sum_i\frac{\mu}{c_i(x(\mu))}\nabla c_i(x(\mu))=0.\tag{17.31} \]

定义乘子估计

\[ \lambda_i(\mu)=\frac{\mu}{c_i(x(\mu))},\tag{17.32} \]

则 \(\nabla f-\sum\lambda_i(\mu)\nabla c_i=0\),即 Lagrange 函数的驻点条件成立;\(c_i>0\)、\(\lambda_i>0\) 也严格成立。KKT 条件中唯一不满足的是互补条件,取而代之的是

\[ \lambda_i(\mu)\,c_i(x(\mu))=\mu.\tag{17.36} \]

所以 \(\mu\downarrow0\) 时 \((x(\mu),\lambda(\mu))\) 越来越接近满足 KKT——这正是第 14 章中心路径的条件 \(x_is_i=\tau\)。

金融直觉:精确的互补条件 \(\lambda_ic_i=0\) 是一个"二选一":要么约束有余量(\(c_i>0\)),影子价格为零;要么约束绑定(\(c_i=0\)),影子价格可以为正。比如一个仓位上限:没顶到上限时,放宽上限毫无价值;顶到了,放宽才值钱。难点在于事先不知道哪些约束会绑定,这是组合性的"猜"。

障碍法把这个硬开关换成 \(\lambda_ic_i=\mu\):每个约束都留一点余量、也都带一点正的影子价格,两者乘积固定为 \(\mu\)。余量越小的约束,估计出的影子价格越大。\(\mu\) 逐步降到零,"谁绑定、谁不绑定"自然浮现,无需事先猜测。这是内点法相对有效集法(第 16a 章)的根本优势。

白话解释:定理 17.4 里的"严格互补"指每个有效约束的乘子都严格为正(不会出现 \(c_i=0\) 且 \(\lambda_i=0\) 的模糊情况);"二阶充分条件"指在约束允许的方向上 Lagrange Hessian 正定,保证 \(x^*\) 是严格局部极小。这几个条件合起来,保证中心路径在 \(x^*\) 附近是一条光滑曲线,可以沿它外推(17.29–17.30)。

定理 17.4:设 \(x^*\) 是局部解,LICQ、严格互补、二阶充分条件成立,\(\mathcal F^o\) 非空。则对充分小的 \(\mu\),存在唯一连续可微的 \(x(\mu)\)(\(x^*\) 邻域内 \(P\) 的局部极小),\(x(\mu)\to x^*\),\(\lambda(\mu)\to\lambda^*\),且 \(\nabla^2_{xx}P\) 正定。轨迹 \(\mathcal C_p=\{x(\mu):\mu>0\}\) 称为原始中心路径。

17.3.5 等式约束与原始–对偶方法

等式约束不能拆成两个不等式(那样 \(\mathcal F^o\) 为空)。办法是对等式加二次罚项:

\[ B(x;\mu)=f(x)-\mu\sum_{i\in\mathcal I}\log c_i(x)+\frac1{2\mu}\sum_{i\in\mathcal E}c_i^2(x).\tag{17.39} \]

要找对不等式严格可行的起点,可引入松弛 \(c_i(x)-s_i=0\)、\(s_i\ge0\),此时任何 \(s>0\) 都在定义域内。

与原始–对偶内点法的关系。原始–对偶方法把乘子 \(\lambda\) 与 \(x\) 同等对待。加松弛 \(s\) 后,障碍问题的一阶条件可写为

\[ \nabla f(x)-A(x)^T\lambda=0,\quad c(x)-s=0,\quad\lambda_is_i=\mu,\quad(\lambda,s)\ge0.\tag{17.41} \]

对数障碍法先用后两式消去 \(s,\lambda\) 再用牛顿法;原始–对偶法则对三个等式整体用(修正)牛顿法:

\[ \begin{bmatrix}\nabla^2_{xx}\mathcal L(x,\lambda)&-A(x)^T&0\\A(x)&0&-I\\0&S&\Lambda\end{bmatrix} \begin{bmatrix}\Delta x\\\Delta\lambda\\\Delta s\end{bmatrix} =\begin{bmatrix}-\nabla f(x)+A(x)^T\lambda\\-c(x)+s\\-\Lambda Se+\mu e+r_{\lambda s}\end{bmatrix},\tag{17.43} \]

\(r_{\lambda s}\) 是修正项(如 Mehrotra 校正)。它与 (14.11)、(14.20)、(16.53) 结构完全相同——LP、QP 的内点步都是 (17.43) 的特例。非线性情形还需要选价值函数、处理非凸性,这是 IPOPT 等软件的核心内容。


17.4 精确罚函数(简介)

二次罚和对数障碍都不精确(需要 \(\mu\to0\))。第 15 章的 \(\ell_1\) 精确罚函数

\[ \phi_1(x;\mu)=f(x)+\frac1\mu\sum_{\mathcal E}|c_i(x)|+\frac1\mu\sum_{\mathcal I}[c_i(x)]^-\tag{17.44} \]

在 \(1/\mu\) 超过最大乘子绝对值时一次极小化就得到解。但它在 \(c_i(x)=0\) 处不可微,而解恰好落在这些点上;一般的不可微优化方法(如丛方法)效率不高。实践中有效的办法是解一列子问题:把 \(f\) 换成以 Lagrange Hessian 为 Hessian 的二次模型、把 \(c_i\) 线性化——这就是第 18 章的 S\(\ell_1\)QP,最强的约束优化技术之一。


17.5 增广拉格朗日法

17.5.1 从罚函数的系统偏差出发

由定理 17.2,二次罚的近似极小点满足

\[ c_i(x_k)\approx-\mu_k\lambda_i^*,\tag{17.45} \]

即存在一个系统性的约束违反,只有 \(\mu_k\to0\) 才消失。增广拉格朗日法的想法是:在目标里加一个线性项把这个偏差"抵消"掉。**增广拉格朗日函数(augmented Lagrangian)**为

\[ \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),\tag{17.46} \]

即 Lagrange 函数加二次罚项。其梯度

\[ \nabla_x\mathcal L_A=\nabla f(x)-\sum_i\Big[\lambda_i-\frac{c_i(x)}{\mu}\Big]\nabla c_i(x).\tag{17.47} \]

固定 \(\lambda^k\)、\(\mu_k\) 极小化 \(\mathcal L_A\) 得 \(x_k\)。与 KKT 条件 \(\nabla f-\sum\lambda_i^*\nabla c_i=0\) 比较,

\[ \lambda_i^*\approx\lambda_i^k-\frac{c_i(x_k)}{\mu_k},\qquad\text{即}\qquad c_i(x_k)\approx-\mu_k(\lambda_i^*-\lambda_i^k).\tag{17.48} \]

关键:若 \(\lambda^k\) 接近 \(\lambda^*\),约束违反量就远小于 \(\mu_k\),而不是与 \(\mu_k\) 成正比。这自然给出乘子更新公式

\[ \lambda_i^{k+1}=\lambda_i^k-\frac{c_i(x_k)}{\mu_k}.\tag{17.49} \]

推导拆解:(17.48) 的来源和 17.2.2 的推导拆解完全相同,只是对号的系数变了。

  1. \(x_k\) 是 \(\mathcal L_A(\cdot,\lambda^k;\mu_k)\) 的近似极小点,所以 (17.47) 近似为零:\(\nabla f(x_k)-\sum_i\big[\lambda_i^k-c_i(x_k)/\mu_k\big]\nabla c_i(x_k)\approx0\)。
  2. KKT 条件:\(\nabla f(x^*)-\sum_i\lambda_i^*\nabla c_i(x^*)=0\)。
  3. 对号得 \(\lambda_i^*\approx\lambda_i^k-c_i(x_k)/\mu_k\),移项即 \(c_i(x_k)\approx-\mu_k(\lambda_i^*-\lambda_i^k)\)。

第 3 步右端是两个小量相乘:\(\mu_k\) 不必小,只要乘子误差 \(\lambda^*-\lambda^k\) 小,违反量就小。(17.49) 则是把第 3 步的近似等号"当真",用它算出下一轮的乘子。

金融直觉:乘子更新像一个不断调价的拍卖师。约束违反 \(c_i(x_k)\) 就是"供需缺口",\(\lambda\) 是价格;每轮按缺口大小调价(步长 \(1/\mu\)),再让参与者(内层优化)按新价格重新决策。价格调对了,缺口自然为零,用不着把罚款单价(\(1/\mu\))抬到天上。二次罚函数相当于价格永远固定为零、只靠罚款逼人就范,所以必须罚到极重才能逼近约束。

框架 17.3(乘子法,等式约束)

给定 μ0 > 0,τ0 > 0,起点 x0^s 和 λ^0;
for k = 0,1,2,...
    从 xk^s 出发求 L_A(·, λ^k; μk) 的近似极小点 xk,当 ‖∇x L_A(x, λ^k; μk)‖ ≤ τk 时停止;
    若满足最终收敛测试:STOP;
    按 (17.49) 更新乘子得 λ^{k+1};
    选新罚参数 μ_{k+1} ∈ (0, μk);
    令 x_{k+1}^s = xk;
end

无需把 \(\mu\) 减到很小就能收敛,病态轻得多,起点也直接用上一步的极小点。Bertsekas 建议内层容差 \(\tau_k=\min(\epsilon_k,\gamma_k\|c(x_k)\|)\),\(\epsilon_k,\gamma_k\to0\)。

例 17.4:对例 17.1,

\[ \mathcal L_A(x,\lambda;\mu)=x_1+x_2-\lambda(x_1^2+x_2^2-2)+\frac1{2\mu}(x_1^2+x_2^2-2)^2.\tag{17.50} \]

取 \(\mu=1\)、\(\lambda=-0.4\):等高线间距说明条件数与 \(Q(x;1)\) 相近,但极小点 \(x_k\approx(-1.02,-1.02)\),比 \(Q(x;1)\) 的 \((-1.1,-1.1)\) 近得多。乘子项带来了实质改进,而且没有付出病态的代价。

17.5.2 理论:为什么 \(\mu\) 不必趋于零

定理 17.5:设 \(x^*\) 是 (17.1) 的局部解,LICQ 成立,\(\lambda=\lambda^*\) 时二阶充分条件成立。则存在 \(\bar\mu>0\),使对所有 \(\mu\in(0,\bar\mu]\),\(x^*\) 是 \(\mathcal L_A(x,\lambda^*;\mu)\) 的严格局部极小。

证明思路。一阶:由 \(c(x^*)=0\),\(\nabla_x\mathcal L_A(x^*,\lambda^*;\mu)=\nabla_x\mathcal L(x^*,\lambda^*)=0\),与 \(\mu\) 无关。二阶:

\[ \nabla^2_{xx}\mathcal L_A(x^*,\lambda^*;\mu)=\nabla^2_{xx}\mathcal L(x^*,\lambda^*)+\frac1\mu A^TA.\tag{17.60} \]

二阶充分条件只保证 \(\nabla^2\mathcal L\) 在零空间(\(Aw=0\))上正定;罚项 \(\frac1\mu A^TA\) 在零空间上为零,但在其正交补(\(A^T\) 的值空间)上提供 \(O(1/\mu)\) 的正曲率。把任意 \(u\) 分解为 \(u=w+A^Tv\)(\(w\in\mathrm{Null}A\),线性代数基本定理),逐项估界并配方得

\[ u^T\nabla^2_{xx}\mathcal L_Au\ge a\big[\|w\|-(b/a)\|v\|\big]^2+\big(d^2/\mu-c-b^2/a\big)\|v\|^2,\tag{17.62} \]

其中 \(a>0\) 是零空间上的曲率下界,\(b,c\) 是交叉项的界,\(d\) 是 \(AA^T\) 的最小特征值。\(\mu\) 足够小使最后一项系数为正即可。∎

白话解释:这个证明的骨架只有一句话:原 Hessian 在"沿约束面"的方向上已经是正的,罚项负责把"离开约束面"的方向也变成正的。 把任意方向 \(u\) 拆成两块:\(w\)(沿约束面,\(Aw=0\))和 \(A^Tv\)(垂直约束面)。在 \(w\) 上,\(\nabla^2\mathcal L\) 提供曲率 \(\ge a\|w\|^2\);在 \(A^Tv\) 上,\(\nabla^2\mathcal L\) 可能是负的(大小不超过 \(c\|v\|^2\)),但罚项提供 \(\frac1\mu\|AA^Tv\|^2\ge\frac{d^2}{\mu}\|v\|^2\) 的正曲率,\(\mu\) 小到一定程度就压得住。交叉项 \(b\) 用配方吸收。(17.62) 的第一项是平方、非负;第二项系数在 \(\mu<d^2/(c+b^2/a)\) 时为正,这就是 \(\bar\mu\)。

注意与二次罚的区别:二次罚要 \(\mu\to0\) 是为了让约束违反消失;这里 \(\mu\le\bar\mu\) 只是为了让 Hessian 正定,一个固定的门槛就够了。

含义:**只要 \(\lambda\) 是 \(\lambda^*\) 的合理估计,即使 \(\mu\) 不很小,极小化 \(\mathcal L_A\) 也能得到 \(x^*\) 的好估计。**第 18 章 SQP 中用增广拉格朗日 Hessian 修正非凸性,也基于这个事实。

定理 17.6(Bertsekas):在定理 17.5 的假设下,存在 \(\delta,\epsilon,M>0\),使对满足 \(\|\lambda^k-\lambda^*\|\le\delta/\mu_k\)、\(\mu_k\le\bar\mu\) 的 \(\lambda^k,\mu_k\):

(a) \(\min_x\mathcal L_A(x,\lambda^k;\mu_k)\) 在 \(\|x-x^*\|\le\epsilon\) 内有唯一解 \(x_k\),且 \(\|x_k-x^*\|\le M\mu_k\|\lambda^k-\lambda^*\|\);

(b) 更新后的乘子满足 \(\|\lambda^{k+1}-\lambda^*\|\le M\mu_k\|\lambda^k-\lambda^*\|\);

(c) \(\nabla^2_{xx}\mathcal L_A(x_k,\lambda^k)\) 正定,\(\nabla c_i(x_k)\) 线性无关。

即乘子误差以因子 \(M\mu_k\) 线性收缩:\(\mu_k\) 固定(只要足够小)就线性收敛,\(\mu_k\) 越小收敛越快。17.7 节的代码在 \(\mu=1\) 固定的情况下观察到每步误差缩小约 9 倍。

17.5.3 不等式约束

方法一:引入松弛 \(c_i(x)-s_i=0\)、\(s_i\ge0\),得到"等式 + 界约束"问题,由 LANCELOT 显式处理界约束(下一小节)。

方法二:在子问题中显式消去松弛(设 \(\mathcal E=\emptyset\))。子问题

\[ \min_{x,s}\ f(x)-\sum_i\lambda_i^k(c_i(x)-s_i)+\frac1{2\mu_k}\sum_i(c_i(x)-s_i)^2\quad\text{s.t.}\quad s_i\ge0\tag{17.52} \]

关于每个 \(s_i\) 是凸二次函数,极小点为

\[ s_i=\max(c_i(x)-\mu\lambda_i^k,\ 0).\tag{17.54} \]

代回得到逐约束的函数

\[ \psi(t,\sigma;\mu)=\begin{cases}-\sigma t+\dfrac1{2\mu}t^2,&t-\mu\sigma\le0\\[2mm]-\dfrac\mu2\sigma^2,&\text{否则},\end{cases}\tag{17.56} \]

推导拆解:固定 \(x\),记 \(t=c_i(x)\)、\(\sigma=\lambda_i^k\),(17.52) 中只和 \(s_i\) 有关的部分是 \(g(s)=\sigma s+\frac1{2\mu}(t-s)^2\)。

  1. 对 \(s\) 求导:\(g'(s)=\sigma-\frac1\mu(t-s)\),令其为零得 \(s=t-\mu\sigma\)。\(g\) 是开口向上的抛物线,所以这是无约束极小点。
  2. 加上 \(s\ge0\):若 \(t-\mu\sigma>0\),取 \(s=t-\mu\sigma\);否则抛物线在 \(s\ge0\) 上单调增,取 \(s=0\)。这就是 (17.54)。
  3. 代回 \(-\sigma(t-s)+\frac1{2\mu}(t-s)^2\):\(s=0\) 时得 \(-\sigma t+\frac{t^2}{2\mu}\);\(s=t-\mu\sigma\) 时 \(t-s=\mu\sigma\),得 \(-\mu\sigma^2+\frac{\mu\sigma^2}2=-\frac\mu2\sigma^2\)。即 (17.56)。

含义:约束余量足够大(\(t>\mu\sigma\))时,\(\psi\) 是与 \(x\) 无关的常数,这个约束对子问题不起作用;余量小或违反时,它按等式约束的增广拉格朗日项起作用。(17.58) 的 \(\max(\cdot,0)\) 则保证不等式乘子非负。

子问题变成无约束的 \(\min_x f(x)+\sum_{i\in\mathcal I}\psi(c_i(x),\lambda_i^k;\mu_k)\),乘子更新为

\[ \lambda_i^{k+1}=\max\Big(\lambda_i^k-\frac{c_i(x_k)}{\mu_k},\ 0\Big).\tag{17.58} \]

\(\psi\) 关于 \(x\) 一阶连续可微,但在 \(c_i(x)=\mu\lambda_i\) 处二阶导不连续。严格互补成立时迭代点通常远离这个区域(有效约束 \(c_i\approx0\) 而 \(\mu\lambda_i\) 明显为正;非有效约束 \(c_i>0\) 而 \(\lambda_i\approx0\))。

17.5.4 LANCELOT 的实用乘子法

LANCELOT(Conn–Gould–Toint)用松弛把不等式变成等式,求解

\[ \min_x f(x)\quad\text{s.t.}\quad c_i(x)=0\ (i=1,\dots,m),\quad l\le x\le u,\tag{17.64} \]

每步子问题是一个界约束问题:

\[ \min_x\ \mathcal L_A(x,\lambda;\mu)\quad\text{s.t.}\quad l\le x\le u,\tag{17.66} \]

用第 16b 章 16.8 节的梯度投影法(配合二次模型和精确或拟牛顿 Hessian)求解。一阶条件写成投影梯度为零:\(P_{[l,u]}\nabla\mathcal L_A=0\)(17.67)。

算法 17.4 的逻辑(省略具体常数):

  • 近似解子问题,使投影梯度 \(\le\omega_k\);
  • 若约束违反 \(\|c(x_k)\|\le\eta_k\)(足够小):认为当前 \(\mu\) 能维持近可行,更新乘子、不减 \(\mu\);
  • 否则:减小 \(\mu\)(\(\mu_{k+1}=\tau\mu_k\)),乘子不变,更重视降低约束违反;
  • 两种情况都收紧 \(\omega_k,\eta_k\),使后续子问题越解越精确。

这种"能更新乘子就不加罚"的策略让 \(\mu\) 通常停在一个适中的值,避免了病态。


17.6 序列线性约束方法(简介)

序列线性约束(SLC)方法,又称约化拉格朗日法,每步在线性化约束下极小化一个非线性目标(SQP 则是二次目标):

\[ \min_x F_k(x)\quad\text{s.t.}\quad\nabla c_i(x_k)^T(x-x_k)+c_i(x_k)=0,\ i\in\mathcal E.\tag{17.69} \]

最流行的选择是增广拉格朗日形式

\[ F_k(x)=f(x)-\sum_i\lambda_i^k\bar c_i^k(x)+\frac1{2\mu}\sum_i[\bar c_i^k(x)]^2,\tag{17.72} \]

其中 \(\bar c_i^k(x)=c_i(x)-c_i(x_k)-\nabla c_i(x_k)^T(x-x_k)\) 是约束与其线性化之差。在满足线性化约束 (17.69) 的点上线性化部分为零,\(\bar c_i^k(x)=c_i(x)\),所以 (17.72) 与增广拉格朗日 (17.46) 一致;用 \(\bar c_i^k\) 代替 \(c_i\) 只是去掉了子问题约束已经保证为零的部分。著名的 MINOS 软件基于此。特点是外迭代很少,但每次外迭代要解一个非线性子问题、函数求值多;普遍认为 SQP 更优,但 MINOS 在大规模、多数约束线性的问题上仍很有效。在量化中直接用到的场景较少。


17.7 量化实战

17.7.1 在哪里用

  • 软约束建模:17.2.3 节的提示——罚系数越大越病态,违反量约为 \(\lambda^*/\rho\)。
  • ADMM:把问题拆成 \(\min f(x)+g(z)\) s.t. \(Ax+Bz=c\),对增广拉格朗日交替极小化 \(x\)、\(z\) 再更新乘子,就是 ADMM。OSQP(大规模 QP)、Lasso / 组 Lasso 因子选择、分布式组合优化都用它。定理 17.5–17.6 解释了为什么它的罚参数可以固定、乘子线性收敛。
  • 风险平价:资产 \(i\) 的风险贡献 \(RC_i=w_i(\Sigma w)_i\)。要求 \(RC_i/\sum_jRC_j=b_i\) 的问题直接写是非凸的,但它有一个障碍型的凸化(第 06 章 6.6 节和第 11 章 11.6.2 节已分别用 Newton–CG 和非线性方程组求解过同一问题,本章补上"它为什么是障碍函数"这一视角):
\[ \min_{y>0}\ \tfrac12y^T\Sigma y-\sum_ib_i\log y_i, \]

驻点条件 \((\Sigma y)_i=b_i/y_i\) 即 \(y_i(\Sigma y)_i=b_i\)——与 (17.36) \(\lambda_ic_i=\mu\) 同构。目标严格凸,有唯一解,归一化 \(w=y/\mathbf 1^Ty\) 后风险贡献比例恰为 \(b\)。

推导拆解:为什么归一化后风险贡献比例恰为 \(b\)?

  1. 梯度:\(\nabla\big(\tfrac12y^T\Sigma y\big)=\Sigma y\),\(\nabla\big(-\sum b_i\log y_i\big)\) 的第 \(i\) 个分量为 \(-b_i/y_i\)。令梯度为零得 \(y_i(\Sigma y)_i=b_i\),即 \(RC_i(y)=b_i\)。
  2. 对 \(i\) 求和:\(\sum_iy_i(\Sigma y)_i=y^T\Sigma y=\sum_ib_i=1\),所以 \(RC_i(y)/\sum_jRC_j(y)=b_i\)。
  3. 风险贡献比例对整体缩放不变:\(w=y/s\)(\(s=\mathbf 1^Ty>0\))时每个 \(RC_i\) 都乘 \(1/s^2\),比例不变。
  4. 严格凸:Hessian \(\Sigma+\mathrm{diag}(b_i/y_i^2)\) 是正定 + 正定,所以驻点唯一,且就是全局极小。

这里对数项起两重作用:作为障碍,保证 \(y>0\)(多头组合);作为"预算",其系数 \(b_i\) 直接成为风险贡献目标。

  • 期权模型校准:SVI、局部波动率校准中的无套利约束和参数界,可用"增广拉格朗日 + 界约束梯度投影"(LANCELOT 式)处理。

17.7.2 代码一:二次罚、增广拉格朗日、对数障碍的数值行为

import numpy as np
from scipy.optimize import minimize

# 例 17.1:min x1+x2 s.t. x1^2+x2^2-2=0;x*=(-1,-1),λ*=-0.5
f  = lambda x: x[0] + x[1]
c  = lambda x: x[0]**2 + x[1]**2 - 2
gf = np.array([1.0, 1.0])
gc = lambda x: 2 * x
def hessQ(x, mu, lam=0.0):      # 增广拉格朗日 Hessian;lam=0 即二次罚 (17.15)
    return (c(x) / mu - lam) * 2 * np.eye(2) + np.outer(gc(x), gc(x)) / mu

print("== 二次罚函数 Q(x;μ) ==")
x = np.array([-1.5, -1.0])
for mu in [1, 0.1, 0.01, 0.001, 1e-4]:
    Q  = lambda x: f(x) + c(x)**2 / (2 * mu)
    gQ = lambda x: gf + c(x) / mu * gc(x)
    x = minimize(Q, x, jac=gQ, method="BFGS", options={"gtol": 1e-10}).x   # 用上一解热启动
    print(f"μ={mu:7.0e}  x(μ)=({x[0]:.5f},{x[1]:.5f})  c(x)={c(x):+.2e}  "
          f"-c/μ={-c(x)/mu:.5f}  cond(∇²Q)={np.linalg.cond(hessQ(x, mu)):.1e}")

print("== 增广拉格朗日(乘子法),固定 μ=1 ==")
mu, lam, x = 1.0, -0.4, np.array([-1.5, -1.0])
for k in range(8):
    LA  = lambda x: f(x) - lam * c(x) + c(x)**2 / (2 * mu)
    gLA = lambda x: gf - (lam - c(x) / mu) * gc(x)
    x = minimize(LA, x, jac=gLA, method="BFGS", options={"gtol": 1e-12}).x
    print(f"k={k}  λ^k={lam:+.7f}  x_k=({x[0]:.6f},{x[1]:.6f})  |λ^k-λ*|={abs(lam+0.5):.1e}  "
          f"cond={np.linalg.cond(hessQ(x, mu, lam)):.1f}")
    lam = lam - c(x) / mu                                   # (17.49)

print("== 对数障碍:例 17.3 ==")
x = np.array([0.5, 0.5])
for mu in [1, 0.1, 0.01, 0.001, 1e-4]:
    P  = lambda x: ((x[0]+0.5)**2 + (x[1]-0.5)**2
                    - mu * np.sum(np.log(x) + np.log(1 - x))) if np.all((x > 0) & (x < 1)) else np.inf
    gP = lambda x: np.array([2*(x[0]+0.5), 2*(x[1]-0.5)]) - mu * (1/x - 1/(1-x))
    hP = lambda x: 2*np.eye(2) + mu * np.diag(1/x**2 + 1/(1-x)**2)
    for _ in range(50):                                     # 带回溯的牛顿法,保持严格可行
        p = -np.linalg.solve(hP(x), gP(x)); t = 1.0
        while P(x + t * p) > P(x) + 1e-4 * t * gP(x) @ p: t *= 0.5
        x = x + t * p
        if np.linalg.norm(gP(x)) < 1e-12: break
    print(f"μ={mu:7.0e}  x(μ)=({x[0]:.6f},{x[1]:.6f})  λ1(μ)=μ/x1={mu/x[0]:.5f}  "
          f"cond(∇²P)={np.linalg.cond(hP(x)):.1e}")

输出:

== 二次罚函数 Q(x;μ) ==
μ=  1e+00  x(μ)=(-1.10716,-1.10716)  c(x)=+4.52e-01  -c/μ=-0.45161  cond(∇²Q)=1.2e+01
μ=  1e-01  x(μ)=(-1.01227,-1.01227)  c(x)=+4.94e-02  -c/μ=-0.49394  cond(∇²Q)=8.4e+01
μ=  1e-02  x(μ)=(-1.00125,-1.00125)  c(x)=+4.99e-03  -c/μ=-0.49938  cond(∇²Q)=8.0e+02
μ=  1e-03  x(μ)=(-1.00012,-1.00012)  c(x)=+5.00e-04  -c/μ=-0.49994  cond(∇²Q)=8.0e+03
μ=  1e-04  x(μ)=(-1.00001,-1.00001)  c(x)=+5.00e-05  -c/μ=-0.49999  cond(∇²Q)=8.0e+04
== 增广拉格朗日(乘子法),固定 μ=1 ==
k=0  λ^k=-0.4000000  x_k=(-1.022059,-1.022059)  |λ^k-λ*|=1.0e-01  cond=9.5
k=1  λ^k=-0.4892086  x_k=(-1.002396,-1.002396)  |λ^k-λ*|=1.1e-02  cond=9.1
k=2  λ^k=-0.4988048  x_k=(-1.000266,-1.000266)  |λ^k-λ*|=1.2e-03  cond=9.0
k=3  λ^k=-0.4998672  x_k=(-1.000030,-1.000030)  |λ^k-λ*|=1.3e-04  cond=9.0
k=4  λ^k=-0.4999852  x_k=(-1.000003,-1.000003)  |λ^k-λ*|=1.5e-05  cond=9.0
k=5  λ^k=-0.4999984  x_k=(-1.000000,-1.000000)  |λ^k-λ*|=1.6e-06  cond=9.0
k=6  λ^k=-0.4999998  x_k=(-1.000000,-1.000000)  |λ^k-λ*|=1.8e-07  cond=9.0
k=7  λ^k=-0.5000000  x_k=(-1.000000,-1.000000)  |λ^k-λ*|=2.0e-08  cond=9.0
== 对数障碍:例 17.3 ==
μ=  1e+00  x(μ)=(0.321037,0.500000)  λ1(μ)=μ/x1=3.11491  cond(∇²P)=1.4e+00
μ=  1e-01  x(μ)=(0.078958,0.500000)  λ1(μ)=μ/x1=1.26649  cond(∇²P)=6.5e+00
μ=  1e-02  x(μ)=(0.009713,0.500000)  λ1(μ)=μ/x1=1.02952  cond(∇²P)=5.2e+01
μ=  1e-03  x(μ)=(0.000997,0.500000)  λ1(μ)=μ/x1=1.00300  cond(∇²P)=5.0e+02
μ=  1e-04  x(μ)=(0.000100,0.500000)  λ1(μ)=μ/x1=1.00030  cond(∇²P)=5.0e+03

三组结果把本章理论逐条兑现:

  1. 二次罚:\(\mu=1\) 时 \(x(\mu)\approx(-1.107,-1.107)\),与原书图 17.1 的 \((-1.1,-1.1)\) 一致。约束违反 \(c(x)\) 与 \(\mu\) 严格成比例(\(c\approx0.5\mu=-\lambda^*\mu\),即 (17.45)),\(-c/\mu\) 收敛到 \(\lambda^*=-0.5\)(定理 17.2),而 Hessian 条件数按 \(1/\mu\) 增长到 \(8\times10^4\)((17.16))。
  2. 增广拉格朗日:\(\mu=1\)、\(\lambda^0=-0.4\) 时第一步就得到 \((-1.022,-1.022)\),与原书例 17.4 的 \((-1.02,-1.02)\) 一致。此后 \(\mu\) 始终固定为 1,乘子误差每步缩小约 9 倍(定理 17.6 的线性收缩),8 步达到 \(10^{-8}\) 精度——条件数始终约为 9,而二次罚要达到同样精度需要 \(\mu\sim10^{-8}\)、条件数 \(\sim10^9\)。
  3. 对数障碍:\(x_1(\mu)\approx\mu\),从内部逼近边界 \(x_1=0\);乘子估计 \(\mu/x_1\to1=\lambda^*\)(\(\nabla f(x^*)=(1,0)=\lambda_1^*e_1\));条件数同样按 \(1/\mu\) 增长((17.27))。带回溯、保持严格可行的牛顿法在每个 \(\mu\) 上都很快收敛,这是框架 17.2 推荐的做法。

17.7.3 代码二:用障碍函数求风险预算组合

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(5)
n = 8                                             # 8 个资产类别
vol = np.array([0.16, 0.18, 0.22, 0.05, 0.07, 0.12, 0.20, 0.03])
C = 0.3 + 0.7 * np.eye(n); C[3:5, 3:5] = [[1, .8], [.8, 1]]; C[0, 3] = C[3, 0] = -0.2
Sigma = np.outer(vol, vol) * C
bud = np.array([0.15, 0.15, 0.10, 0.15, 0.15, 0.10, 0.10, 0.10])   # 风险预算

# 障碍型凸化:min 1/2 y'Σy - Σ b_i log y_i  (无约束、严格凸,y>0 自动保持)
# 驻点条件 (Σy)_i = b_i / y_i  ⇔  y_i (Σy)_i = b_i ,与 (17.36) λ_i c_i = μ 同构
phi  = lambda y: 0.5 * y @ Sigma @ y - bud @ np.log(y) if np.all(y > 0) else np.inf
grad = lambda y: Sigma @ y - bud / y
hess = lambda y: Sigma + np.diag(bud / y**2)
y = np.full(n, 1.0)
for k in range(30):                               # 牛顿法 + 保持 y>0 的回溯
    g = grad(y)
    if np.linalg.norm(g) < 1e-13: break
    p = -np.linalg.solve(hess(y), g); t = 1.0
    while phi(y + t * p) > phi(y) + 1e-4 * t * g @ p: t *= 0.5
    y = y + t * p
w = y / y.sum()
rc = w * (Sigma @ w) / (w @ Sigma @ w)            # 风险贡献占比
print("牛顿迭代次数:", k)
print("权重     :", w.round(4))
print("风险贡献 :", rc.round(4))
print("目标预算 :", bud)
print(f"组合波动 : {np.sqrt(w @ Sigma @ w):.2%}")

# 对照:直接用 SLSQP 求解"风险贡献 = 预算"的非凸最小二乘形式
def rc_err(w):
    s = Sigma @ w; tot = w @ s
    return np.sum((w * s / tot - bud) ** 2)
r = minimize(rc_err, np.full(n, 1 / n), method="SLSQP", bounds=[(1e-6, 1)] * n,
             constraints=[{"type": "eq", "fun": lambda w: w.sum() - 1}],
             options={"ftol": 1e-16, "maxiter": 500})
print(f"SLSQP 直接法: 迭代 {r.nit},与障碍法权重最大差 {np.abs(r.x - w).max():.1e}")

输出:

牛顿迭代次数: 7
权重     : [0.0882 0.067  0.0393 0.2521 0.1493 0.0721 0.0433 0.2886]
风险贡献 : [0.15 0.15 0.1  0.15 0.15 0.1  0.1  0.1 ]
目标预算 : [0.15 0.15 0.1  0.15 0.15 0.1  0.1  0.1 ]
组合波动 : 5.20%
SLSQP 直接法: 迭代 20,与障碍法权重最大差 5.7e-08

障碍型凸化只用 7 步牛顿迭代就得到风险贡献与预算完全一致的组合;直接对"风险贡献 = 预算"做非凸最小二乘也能得到同一解,但需要更多迭代,而且在资产多、预算极不均衡时容易卡在局部解或边界上。凸化的理论保证来自本章:对数项使 \(y\) 自动保持为正(障碍性质),驻点条件是"扰动互补"形式,严格凸保证唯一性。低波动资产(第 4、8 类,波动 5% 和 3%)分到了最大的权重——这是风险平价"加杠杆低波资产"特征的来源。


本章小结

罚、障碍和增广拉格朗日方法都用一列无约束(或界约束)子问题逼近约束问题。二次罚函数 \(f+\frac1{2\mu}\|c\|^2\) 的极小点在 \(\mu\to0\) 时收敛到 KKT 点,\(-c_i/\mu\) 估计乘子,但约束违反量与 \(\mu\) 成正比、Hessian 条件数按 \(1/\mu\) 增长。对数障碍 \(f-\mu\sum\log c_i\) 让迭代点保持严格可行,\(\mu/c_i\) 估计乘子,\(\lambda_ic_i=\mu\) 是扰动互补,极小点轨迹就是中心路径;凸问题全局收敛,良态局部解附近存在光滑路径;同样病态,需要牛顿法加外推起点。把扰动 KKT 系统直接用牛顿法求解就是原始–对偶内点法,(17.43) 统一了 LP、QP、NLP 的内点步。\(\ell_1\) 精确罚只需一次极小化但不光滑,实用中通过 S\(\ell_1\)QP 实现。增广拉格朗日 \(f-\lambda^Tc+\frac1{2\mu}\|c\|^2\) 配合乘子更新 \(\lambda\leftarrow\lambda-c/\mu\),在 \(\lambda=\lambda^*\) 时对一切小于阈值的 \(\mu\) 精确,乘子误差以 \(M\mu\) 线性收缩,所以 \(\mu\) 不必趋于零、没有严重病态;不等式通过 \(\psi\) 函数或松弛 + 界约束(LANCELOT)处理。SLC/MINOS 在线性化约束下极小化增广拉格朗日函数。

概念 公式 / 要点
二次罚 \(Q=f+\frac1{2\mu}\sum c_i^2\);\(c_i(x_k)\approx-\mu_k\lambda_i^*\)
罚的 Hessian \(\nabla^2Q\approx\nabla^2\mathcal L+\frac1\mu A^TA\),条件数 \(O(1/\mu)\)
稳定牛顿方程 (17.18):增广 \([\cdot\ A^T;\ A\ -\mu I]\)
对数障碍 \(P=f-\mu\sum\log c_i\);\(\lambda_i(\mu)=\mu/c_i\);\(\lambda_ic_i=\mu\)
原始–对偶步 (17.43),LP/QP 为其特例
增广拉格朗日 \(\mathcal L_A=f-\lambda^Tc+\frac1{2\mu}|c|^2\)
乘子更新 \(\lambda^{k+1}=\lambda^k-c(x_k)/\mu_k\);不等式取 \(\max(\cdot,0)\)
定理 17.5 \(\nabla^2\mathcal L_A=\nabla^2\mathcal L+\frac1\mu A^TA\) 在 \(\mu\le\bar\mu\) 时正定
定理 17.6 \(|\lambda^{k+1}-\lambda^*|\le M\mu_k|\lambda^k-\lambda^*|\)
\(\psi\) 函数 \(\psi(t,\sigma;\mu)=-\sigma t+t^2/(2\mu)\)(\(t\le\mu\sigma\)),否则 \(-\mu\sigma^2/2\)
风险平价凸化 \(\min\frac12y^T\Sigma y-\sum b_i\log y_i\),\(w=y/\mathbf 1^Ty\)

练习

基础

  1. 证明 \(\min(0,z)^2\) 在 \(z=0\) 处一阶可导但二阶导不连续,并说明这对 (17.5) 的含义。
  2. 对 \(\min x\) s.t. \(x\ge0\),写出二次罚函数 \(\frac1{2\mu}\min(0,x)^2+x\) 的极小点,验证 \(-c/\mu\to\lambda^*=1\)。
  3. 证明 \(\min\frac1{1+x^2}\) s.t. \(x\ge1\) 的对数障碍函数对任意 \(\mu>0\) 无下界。
  4. 证明 (17.18) 与 \(\nabla^2_{xx}Q\,p=-\nabla_xQ\) 给出相同的 \(p\)。
  5. 对 \(\psi(t,\sigma;\mu)\) 证明它在 \(t=\mu\sigma\) 处连续可微,写出两侧的二阶导。

进阶

  1. 对 \(\min x\) s.t. \(x^2\ge0\),\(x+1\ge0\)(解 \(x^*=-1\)),证明对数障碍函数存在收敛到 \(0\)(非解)的局部极小序列。
  2. 从 \(P(\cdot;\mu_k)\) 的精确极小点出发对 \(P(\cdot;\mu_{k+1})\) 做一步牛顿,与 (17.29)–(17.30) 的切向预测比较。
  3. 把 (17.41) 推广到同时含等式与不等式约束的问题,写出对应的原始–对偶牛顿步。
  4. 用增广拉格朗日法(框架 17.3 + 不等式的 \(\psi\) 函数)求解第 16b 章中的长仓最小方差组合(预算等式 + \(w\ge0\)),内层用 BFGS,比较 \(\mu\) 固定为 \(1,0.1,0.01\) 时的外迭代次数。
  5. 实现 ADMM 求解 \(\min\frac12w^T\Sigma w-\gamma\alpha^Tw\) s.t. \(\mathbf 1^Tw=1\),\(0\le w\le u\):令 \(z=w\),\(x\)-步解带预算约束的二次问题,\(z\)-步是截断投影,再更新乘子。与 SLSQP 的解比较,并观察罚参数对收敛速度的影响。

原书推荐习题:17.1、17.8(罚函数 / 增广拉格朗日的光滑性)、17.2、17.3(障碍法的失败模式)、17.5(中心路径外推与牛顿预测步)、17.6(含等式约束的原始–对偶系统)、17.7(投影梯度一阶条件)。原书第 17 章共 17.1–17.8 题。


原书对照

本章内容 原书位置 PDF 页码
17.1 核心思路 第 17 章引言 PDF p.506–508
17.2 二次罚函数、框架 17.1、定理 17.1–17.2、病态与 (17.18) §17.1 The Quadratic Penalty Method,式 (17.1)–(17.18) PDF p.508–516
17.3 对数障碍、例 17.2–17.3、定理 17.3–17.4、原始–对偶关系 §17.2 The Logarithmic Barrier Method,式 (17.19)–(17.43) PDF p.516–528
17.4 精确罚函数 §17.3 Exact Penalty Functions,式 (17.44) PDF p.528–529
17.5 增广拉格朗日、例 17.4、定理 17.5–17.6、\(\psi\) 函数、LANCELOT §17.4 Augmented Lagrangian Method,式 (17.45)–(17.68),算法 17.4 PDF p.529–540
17.6 序列线性约束方法 §17.5 Sequential Linearly Constrained Methods,式 (17.69)–(17.72) PDF p.540–541
注释与参考 Notes and References PDF p.541–542
习题 17.1–17.8 Exercises PDF p.542–543

说明:本章结构与原书第 1 版一致(对数障碍法在第 17 章;第 2 版移到第 19 章)。定理 17.5 证明中"代数基本定理"应为"线性代数基本定理"。原书正文页码 = PDF 页码减 18(已用 PDF 页眉核对,第 17 章至附录均如此)。