第 16b 章二次规划：梯度投影法、内点法、对偶与组合优化实战

本篇接第 16a 章（16a_二次规划_KKT系统与有效集法.md）。上篇的有效集方法每步只改一个约束，在"解处有上百个约束有效"的大问题上会很慢。本篇介绍两类能处理大规模问题的方法——梯度投影法和内点法，再讲 QP 对偶，最后用 scipy 完整求解一个带换手、行业、个股上限约束的组合优化问题。

学习目标

读完本章下篇，你应当能够：

说明梯度投影法为什么能一次大幅改变有效集；会计算界约束 QP 的 Cauchy 点（沿投影梯度分段路径的第一个局部极小），并理解子空间极小化步。
把第 14 章的原始–对偶内点法推广到凸 QP：写出 KKT 条件、中心路径、牛顿步方程 (16.53)，以及增广系统与正规方程两种消元形式；知道 QP 中原始、对偶必须用同一步长。
写出严格凸 QP 的对偶问题，理解它为什么是一个界约束 QP。
比较有效集法、梯度投影法、内点法的适用场景，并能为具体的组合优化问题选择求解器。
用 scipy（SLSQP、trust-constr、linprog、L-BFGS-B）和自写内点法求解含预算、换手拆分、个股上限、行业偏离、换手预算的均值–方差组合，读出并验证影子价格。

读前导读

这一章在解决什么问题

上篇的有效集方法一次只能把一只股票"钉上"或"放开"。几十只股票时没问题；几千只股票、解处有几百只顶在限额上时，就要几百轮迭代。本篇讲两种能"成批调整"的方法，再讲 QP 的对偶，最后把一个真实感很强的组合问题从头到尾解一遍。

梯度投影法的思路就是组合经理常做的"先按信号调、再截断到限额"：沿最有利的方向同时调整所有权重，碰到上下限的股票就停在限额上，其他股票继续调。一步之内可以有上百只股票同时撞到限额，所以特别适合"只有个股上下限"的多空组合。内点法则是第 14 章 LP 内点法的直接推广，迭代次数几十次，几乎不随股票数增长。

本篇的重头戏是 16.11 节的完整实例：带预算、换手拆分、个股上限、行业偏离、换手预算的均值–方差组合。你会看到每条约束的乘子怎样变成可以拿去和投资经理、风控讨论的数字。例如换手预算的乘子回答"多给 1% 换手额度，组合效用能提高多少"，这和第 12 章资本预算中"多给 1 元预算能多赚多少 NPV"是同一个问题。16.11.2 节后的讲解框还会推出一个实务上很有用的结论：交易成本和换手预算会在每只股票周围形成一个"不交易区间"。

需要先想起来的数学

1. 投影（截断）。 把一个点"投影"到盒子 \([l,u]\) 上，就是逐分量截断：低于下限取下限，高于上限取上限，中间的不动。例：上限 2%、下限 −2% 时，\((3\%,-1\%,-5\%)\) 投影为 \((2\%,-1\%,-2\%)\)。Excel 里就是 MAX(l, MIN(u, x))。

2. 一元二次函数的极小点。 \(f(t)=a+bt+\tfrac12ct^2\)（\(c>0\)）的极小点在 \(t^*=-b/c\)。16.8.2 节沿每一段路径做的就是这件事。见第 00 册第 02 章导数与泰勒展开。

3. 共轭梯度法（CG）。 解 \(Gx=-d\)（\(G\) 对称正定）的迭代方法，只需要做矩阵乘向量，最多 \(n\) 步得到精确解，通常远少于 \(n\) 步就足够好。本章只把它当作"求解自由变量子问题的工具"，不需要掌握细节（本册前面的章节有专门介绍）。

4. 对偶与 KKT。 第 13 章的 LP 对偶（对偶变量 = 约束的价格）和第 12 章的 KKT 条件、影子价格在本章全部用到；第 14 章的中心路径、\(\mu\)、Mehrotra 预测–校正在 16.9 节直接照搬。建议先确认这三处已经读懂。见第 00 册第 05 章多元微积分与优化。

怎么读这一章

核心必读是 16.8.1（梯度投影的动机）、16.9.4（有效集法与内点法的对比表）、16.10（QP 对偶的结论与量化提示）和 16.11 全部（组合优化实例，尤其是 16.11.2 的逐条解读和 16.11.4 的实务清单）。

第一次可以只看结论的是 16.8.2 的 Cauchy 点计算细节、16.8.3 的子空间极小化和 16.9.2–16.9.3 的牛顿步与线性代数（读过第 14 章的话会发现结构完全相同）。如果时间有限，可以直接从 16.11 读起，遇到不懂的概念再往回翻。

16.8 梯度投影法

16.8.1 动机

经典有效集法每步只增删一个约束。若初始点没有有效约束、而解处有 200 个有效约束，至少需要 200 次迭代。**梯度投影法（gradient-projection method）**能在一步之内改变许多约束的状态，在约束很简单——尤其是只有变量上下界——时效率最高。

考虑界约束 QP：

\min_x\ q(x)=\tfrac12x^TGx+x^Td\quad\text{s.t.}\quad l\le x\le u.\tag{16.43}

可行域是一个"盒子"，缺失的界取 \(\pm\infty\)。不要求 \(G\) 正定，凸与非凸都可用。这类问题比看上去常见：下面 16.10 节会看到，严格凸 QP 的对偶就是界约束 QP；第 17 章的增广拉格朗日法（LANCELOT）每步要解的也是界约束子问题。

每次迭代分两阶段：

从当前点沿最速下降方向 \(-g\)（\(g=Gx+d\)）出发，碰到界就把方向"折弯"以保持可行，沿这条分段线性路径找 \(q\) 的第一个局部极小点，称为 Cauchy 点 \(x^c\)（类比原书第 4 章信赖域的 Cauchy 点）。把在 \(x^c\) 处有效的界作为工作集。
固定这些有效的界，在剩余的"自由变量"上进一步降低 \(q\)（子空间极小化）。

16.8.2 Cauchy 点

投影算子把点逐分量截断到盒子里：

P(x,l,u)_i=\begin{cases}l_i,&x_i<l_i\\x_i,&l_i\le x_i\le u_i\\u_i,&x_i>u_i.\end{cases}\tag{16.44}

分段线性路径为 \(x(t)=P(x^0-tg,l,u)\)（16.45）。第 \(i\) 个分量撞到界的时刻是

\bar t_i=\begin{cases}(x_i^0-u_i)/g_i,&g_i<0,\ u_i<+\infty\\(x_i^0-l_i)/g_i,&g_i>0,\ l_i>-\infty\\\infty,&\text{其他}\end{cases}\tag{16.46}

白话解释：投影路径 \(x(t)=P(x^0-tg,l,u)\) 可以这样想象：每只股票的权重按自己的"边际吸引力" \(-g_i\) 匀速变化（吸引力为正就加仓、为负就减仓），\(t\) 是时间。某只股票撞到上限或下限就停在那里，其余股票继续走。(16.46) 的 \(\bar t_i\) 就是第 \(i\) 只股票撞线的时刻：距离（到限额还有多远）除以速度 \(|g_i|\)。把撞线时刻排序，就得到一串"断点"，两个断点之间仍在移动的股票集合不变，路径是一段直线。和有效集法每轮只钉住一只股票相比，这里一条路径上可以先后钉住很多只。

每个分量以速率 \(-g_i\) 移动，碰界后停住。把 \(\bar t_i\) 去重排序为断点 \(0<t_1<t_2<\cdots\)，依次检查每段 \([t_{j-1},t_j]\)。在该段上 \(x(t)=x(t_{j-1})+\Delta t\,p^{j-1}\)，方向

p_i^{j-1}=\begin{cases}-g_i,&t_{j-1}<\bar t_i\\0,&\text{其他}\end{cases}\tag{16.47}

\(q\) 沿该段是 \(\Delta t\) 的一元二次函数：

q(x(t))=f_{j-1}+f'_{j-1}\Delta t+\tfrac12f''_{j-1}(\Delta t)^2,\tag{16.48}

\(f'_{j-1}=d^Tp^{j-1}+x(t_{j-1})^TGp^{j-1}\)，\(f''_{j-1}=(p^{j-1})^TGp^{j-1}\)。规则：

若 \(f'_{j-1}\ge0\)，本段起点就是局部极小，\(x^c=x(t_{j-1})\)；
否则若 \(f''_{j-1}>0\) 且 \(\Delta t^*=-f'_{j-1}/f''_{j-1}<t_j-t_{j-1}\)，极小点在段内 \(t_{j-1}+\Delta t^*\)；
否则进入下一段。

相邻两段的方向通常只差一个分量，系数可以增量更新，所以找 Cauchy 点的代价很低。

推导拆解：(16.48) 的系数怎么来。在一段上 \(x=x(t_{j-1})+\Delta t\,p\)（记 \(p=p^{j-1}\)，\(\bar x=x(t_{j-1})\)），代入 \(q(x)=\tfrac12x^TGx+x^Td\) 并按 \(\Delta t\) 的幂次整理：常数项 \(\tfrac12\bar x^TG\bar x+\bar x^Td=q(\bar x)=f_{j-1}\)；一次项 \(\Delta t\,(\bar x^TGp+d^Tp)\)，括号就是 \(f'_{j-1}\)，也等于 \(\nabla q(\bar x)^Tp\)，即"沿这个方向的斜率"（\(G\) 对称，所以 \(\tfrac12\bar x^TGp+\tfrac12p^TG\bar x=\bar x^TGp\)）；二次项 \(\tfrac12(\Delta t)^2p^TGp\)，所以 \(f''_{j-1}=p^TGp\)，即沿这个方向的曲率。三条规则就是一元二次函数求极小：斜率已经非负说明起点就是最低点；斜率为负且向上弯（\(f''>0\)）时极小点在 \(-f'/f''\)，若它落在本段内就停下；否则说明在本段内一直下降，走到段末进入下一段。

16.8.3 子空间极小化

得到 \(x^c\) 后，在有效界固定的子空间上改进：

\min_x q(x)\quad\text{s.t.}\quad x_i=x_i^c\ (i\in\mathcal A(x^c)),\quad l_i\le x_i\le u_i\ (i\notin\mathcal A(x^c)).\tag{16.49}

它不必（也不宜）精确求解——它几乎和原问题一样难。全局收敛只要求近似解 \(x^+\) 可行且 \(q(x^+)\le q(x^c)\)。常用折中：暂时忽略自由变量的界，从 \(x^c\) 出发对自由变量做 CG，一旦碰到界或遇到负曲率就停止（类似原书算法 4.3 的 Steihaug–CG）。这个子问题的零空间基极其简单——就是自由变量对应的单位向量。

算法 16.2（QP 的梯度投影法）

计算可行初始点 x0；
for k = 0,1,2,...
    若 xk 满足 (16.43) 的 KKT 条件：STOP；
    令 x = xk，求 Cauchy 点 xc；
    求 (16.49) 的近似解 x+，使 q(x+) ≤ q(xc) 且可行；
    x_{k+1} = x+；
end

收敛性质：若严格互补成立（解处有效界的乘子都非零），\(\mathcal A(x^c)\) 最终稳定，不再反复进出；退化时可能不稳定。

一般线性约束：原理上也可用，但投影到 \(\{a_i^Tx\ge b_i\}\) 本身就要解一个凸 QP \(\min\|x-\bar x\|^2\) s.t. \(a_i^Tx\ge b_i\)，通常不划算。（一个重要例外是单纯形 \(\{x\ge0,\ \mathbf 1^Tx=1\}\) 和"盒子 + 一条等式"，它们的投影有 \(O(n\log n)\) 的排序算法——这正是"长仓 + 预算"组合的约束集。）

16.9 凸 QP 的内点法

把第 14 章 LP 的原始–对偶内点法推广到凸 QP，描述简单、相对易实现，在大规模问题上很有效。考虑

\min_x\ q(x)=\tfrac12x^TGx+x^Td\quad\text{s.t.}\quad Ax\ge b,\tag{16.50}

\(G\) 对称半正定，\(A\) 为 \(m\times n\)。等式约束只需简单修改（见 16.11 节的实现）。

16.9.1 KKT 条件与中心路径

KKT 条件：\(Gx-A^T\lambda+d=0\)，\(Ax-b\ge0\)，\((Ax-b)_i\lambda_i=0\)，\(\lambda\ge0\)。引入松弛 \(y=Ax-b\)：

Gx-A^T\lambda+d=0,\qquad Ax-y-b=0,\qquad y_i\lambda_i=0,\qquad (y,\lambda)\ge0.\tag{16.51}

与 LP 的 KKT (14.3) 结构相同。由于目标和可行域都凸，KKT 条件既必要又充分。对偶度量

\mu=\frac1m\sum_{i=1}^my_i\lambda_i=\frac{y^T\lambda}m,\tag{16.52}

中心路径是把互补条件换成 \(y_i\lambda_i=\tau\) 的严格可行点曲线。

16.9.2 牛顿步

朝中心路径上 \(\tau=\sigma\mu\) 的点走牛顿步：

\begin{bmatrix}G&0&-A^T\\A&-I&0\\0&\Lambda&Y\end{bmatrix} \begin{bmatrix}\Delta x\\\Delta y\\\Delta\lambda\end{bmatrix} =\begin{bmatrix}-r_d\\-r_b\\-\Lambda Ye+\sigma\mu e\end{bmatrix},\tag{16.53}

\(r_d=Gx-A^T\lambda+d\)，\(r_b=Ax-y-b\)，\(Y=\mathrm{diag}(y)\)，\(\Lambda=\mathrm{diag}(\lambda)\)。更新 \((x,y,\lambda)+\alpha(\Delta x,\Delta y,\Delta\lambda)\)，\(\alpha\) 保证 \((y,\lambda)>0\)。

Mehrotra 预测–校正同样适用，唯一的例外是：原始变量 \((x,y)\) 与对偶变量 \(\lambda\) 必须用同一个步长。LP 中原始、对偶方程不耦合，可以各取各的步长；QP 中它们通过 \(G\) 耦合在 \(Gx-A^T\lambda+d=0\) 里，不同步长会破坏这个方程的可行性。

推导拆解：为什么必须同一步长。牛顿方程第一行保证 \(G\Delta x-A^T\Delta\lambda=-r_d\)。若 \(x\) 走 \(\alpha_p\)、\(\lambda\) 走 \(\alpha_d\)，新残差是 \(r_d+\alpha_pG\Delta x-\alpha_dA^T\Delta\lambda\)。\(\alpha_p=\alpha_d=\alpha\) 时它等于 \((1-\alpha)r_d\)，按比例缩小，全步时归零；两者不等时，残差里会多出 \((\alpha_p-\alpha_d)G\Delta x\) 这样一项，不再保证下降。LP 中 \(G=0\)，第一行 \(A^T\lambda+s=c\) 只含对偶变量，第二行 \(Ax=b\) 只含原始变量，两边各自按比例缩小，互不干扰，所以可以分别取步长。这里的松弛变量 \(y=Ax-b\) 就是"约束还剩多少额度"，互补条件 \(y_i\lambda_i=0\) 仍是那句话：还有剩余额度的约束，影子价格为零。

16.9.3 线性代数

消去 \(\Delta y=A\Delta x+r_b\)，得到增广系统

\begin{bmatrix}G&-A^T\\A&\Lambda^{-1}Y\end{bmatrix} \begin{bmatrix}\Delta x\\\Delta\lambda\end{bmatrix} =\begin{bmatrix}-r_d\\-r_b+(-y+\sigma\mu\Lambda^{-1}e)\end{bmatrix},\tag{16.54}

可用对称不定分解。再消去 \(\Delta\lambda\)，得到正规方程

\big(G+A^T(Y^{-1}\Lambda)A\big)\Delta x=-r_d+A^T(Y^{-1}\Lambda)\big[-r_b-y+\sigma\mu\Lambda^{-1}e\big],

系数矩阵对称正半定（通常正定），可用（修正）Cholesky。每步 \(y,\lambda\) 变化会改变 \(A^T(Y^{-1}\Lambda)A\) 的数值，所以每步都要重新分解，但非零结构不变，可以重用符号分解。

邻域和势函数也可以照搬：\(\mathcal N_2(\theta)\)、\(\mathcal N_{-\infty}(\gamma)\) 中把 \(x_is_i\) 换成 \(y_i\lambda_i\)；势函数 \(\Phi_\rho=\rho\log y^T\lambda-\sum\log y_i\lambda_i\)，\(\rho>m\)。

16.9.4 有效集法与内点法的对比

	有效集法	内点法
迭代次数	多（至少等于有效约束数的变化量）	少（几十次，几乎与规模无关）
每步代价	低（更新分解）	高（重新分解）
实现难度	高，尤其稀疏分解的更新	较低，可直接用标准稀疏 Cholesky
热启动	很好（昨日工作集）	困难
解的形式	精确落在有效约束上	内部近似解，需要阈值/交叉清理
结构	几次基更新后带状、稀疏结构就丢失	块带状（最优控制、多期组合）可完整利用

16.10 QP 的对偶

\(G\) 正定时，

\min_x\ \tfrac12x^TGx+x^Td\quad\text{s.t.}\quad Ax\ge b\tag{16.55}

的（Wolfe）对偶为

\max_{x,\lambda}\ \tfrac12x^TGx+x^Td-\lambda^T(Ax-b)\quad\text{s.t.}\quad Gx+d-A^T\lambda=0,\ \lambda\ge0.\tag{16.56}

由约束解出 \(x=G^{-1}(A^T\lambda-d)\) 代入，得到只有非负约束的问题

\max_{\lambda\ge0}\ -\tfrac12\lambda^T(AG^{-1}A^T)\lambda+\lambda^T(b+AG^{-1}d)-\tfrac12d^TG^{-1}d.\tag{16.57}

这是一个界约束 QP，可以用梯度投影法求解，通常比经典有效集法更快识别有效集。缺点是目标中的 \(AG^{-1}A^T\) 形式复杂：直接法需要显式形成它；替代做法是分解 \(G\) 后用 CG。

金融直觉：对偶问题可以读成"给约束定价，然后放手让交易员自由优化"。给定一组约束价格 \(\lambda\ge0\)，把"消耗约束额度的代价" \(-\lambda^T(Ax-b)\) 并入目标后，问题就没有约束了，最优组合直接写成 \(x(\lambda)=G^{-1}(A^T\lambda-d)\)，这就是 (16.56) 的等式约束。价格定得太低，自由优化的结果会超限；定得太高，又白白放弃收益。对偶问题 (16.57) 就是在所有非负价格中找"恰好让自由优化的结果守住限额"的那一组，它就是原问题的乘子 \(\lambda^*\)。这与机构内部"风险预算定价"的做法同理：风控不逐笔审批，而是给每单位风险额度定一个内部价格，交易台在扣除这个价格后自由决策；价格定对了，各台加总的结果自然满足总限额。 (16.57) 的推导只用了代入：由 \(Gx=A^T\lambda-d\)，Lagrange 函数中 \(x^T(d-A^T\lambda)=-x^TGx\)，于是 \(\tfrac12x^TGx+x^Td-\lambda^TAx+\lambda^Tb=-\tfrac12x^TGx+\lambda^Tb\)，再把 \(x=G^{-1}(A^T\lambda-d)\) 代入展开即得。

量化提示：对偶变量 \(\lambda\) 的维数是约束个数，原始变量的维数是资产个数。约束少、资产多（如数千只股票、几十条因子/行业约束、无个股界）时，对偶问题小得多；因子模型下 \(G^{-1}=\Sigma^{-1}\) 还可以用 Sherman–Morrison–Woodbury 快速计算（附录 A）。

16.11 量化实战（下）：均值–方差组合的完整求解

16.11.1 模型

给定 40 只股票的 alpha 预测 \(\alpha\)、因子模型协方差 \(\Sigma=BFB^T+D\)、当前持仓 \(w^0\)、基准 \(w_b\)、4 个行业。决策变量 \(x=(w,b,s)\in\mathbb R^{3n}\)，其中 \(b,s\) 是买入量和卖出量：

\begin{aligned} \min_{w,b,s}\quad&\kappa\,w^T\Sigma w-\alpha^Tw+c\,\mathbf 1^T(b+s)\\ \text{s.t.}\quad&\mathbf 1^Tw=1,\qquad w-b+s=w^0 &&\text{（预算、换手拆分：等式）}\\ &0\le w\le u,\quad b\ge0,\ s\ge0 &&\text{（长仓、个股上限）}\\ &|H^T(w-w_b)|\le\delta &&\text{（行业偏离 ±2\%）}\\ &\mathbf 1^T(b+s)\le\tau &&\text{（换手预算）} \end{aligned}

参数：\(\kappa=2\)，单边费率 \(c=20\)bp，\(u=8\%\)，\(\delta=2\%\)。这是一个凸 QP：\(G=\mathrm{diag}(2\kappa\Sigma,0,0)\) 半正定（在 \(b,s\) 方向上是零曲率，由线性成本和非负约束控制）。

几个建模要点：

换手拆分是精确的：成本 \(c>0\) 使最优解不会同时买卖同一只股票（\(b_is_i=0\)），所以 \(\mathbf 1^T(b+s)=\|w-w^0\|_1\)。

推导拆解：为什么 \(b_is_i=0\) 自动成立。假设某个解里同一只股票既买 \(b_i=3\%\) 又卖 \(s_i=1\%\)，净变化 \(w_i-w_i^0=2\%\)。把两者同时减去 \(\min(b_i,s_i)=1\%\)，变成只买 2%、不卖：\(w\) 完全不变，所以风险、alpha、预算、行业约束都不受影响，但交易成本少了 \(2c\times1\%\)，换手也少了 2%（换手预算只会更宽松）。所以任何"既买又卖"的解都能被严格改进，最优解里不可能出现。若 \(c=0\)，这个论证失效，见练习 5。

可行性先于最优性：行业偏离和个股上限要求从 \(w^0\) 出发必须做一定量的交易。换手预算低于这个最小值时问题不可行。代码先用一个 LP（第 13 章）算出最小换手。
求解器：scipy 没有专门的 QP 求解器。SLSQP 是线搜索 SQP（第 18 章），对凸 QP 它每步解的就是原问题本身的线性化——一般几步就收敛；trust-constr 对不等式用障碍法内核（第 17 章）；我们再自写一个 16.9 节的原始–对偶内点法作为对照。实务中常用 OSQP、Clarabel、MOSEK、Gurobi（通过 CVXPY 建模），原理都在本章。

16.11.2 代码一：三种求解器、影子价格与换手预算扫描

import numpy as np, time
from scipy.optimize import minimize, LinearConstraint, Bounds

# ---------------- 1. 模拟数据：因子模型协方差 + 行业 + 当前持仓 ----------------
rng = np.random.default_rng(2024)
n, nf, nind = 40, 3, 4
Bexp = rng.standard_normal((n, nf)) * np.array([0.9, 0.4, 0.3]) + np.array([1.0, 0, 0])
Fcov = np.diag([0.04, 0.01, 0.008])                    # 年化因子协方差
Dvar = rng.uniform(0.02, 0.09, n)                       # 特质方差
Sigma = Bexp @ Fcov @ Bexp.T + np.diag(Dvar)
alpha = 0.04 * rng.standard_normal(n)                   # 年化 alpha 预测
ind = rng.integers(0, nind, n)
H = np.array([(ind == j).astype(float) for j in range(nind)])
w_bench = rng.dirichlet(np.full(n, 2.0))                # 基准权重
w0 = rng.dirichlet(np.full(n, 2.0))                     # 当前持仓
kappa, tcost, ub, ind_dev = 2.0, 0.002, 0.08, 0.02      # 风险厌恶、单边费率、个股上限、行业偏离

# 变量 x = [w; b; s]，w = w0 + b - s，b,s >= 0 分别为买入、卖出
N = 3 * n
G = np.zeros((N, N)); G[:n, :n] = 2 * kappa * Sigma
dvec = np.r_[-alpha, np.full(n, tcost), np.full(n, tcost)]
Aeq = np.vstack([np.r_[np.ones(n), np.zeros(2 * n)],
                 np.hstack([np.eye(n), -np.eye(n), np.eye(n)])])
beq = np.r_[1.0, w0]
def ineq(tau):
    """A_in x >= b_in：w>=0, -w>=-ub, b>=0, s>=0, 行业上下限, 换手上限"""
    I, Z = np.eye(n), np.zeros((n, n))
    Ain = np.vstack([np.hstack([I, Z, Z]), np.hstack([-I, Z, Z]),
                     np.hstack([Z, I, Z]), np.hstack([Z, Z, I]),
                     np.hstack([H, np.zeros((nind, 2 * n))]),
                     np.hstack([-H, np.zeros((nind, 2 * n))]),
                     np.r_[np.zeros(n), -np.ones(2 * n)][None, :]])
    bin_ = np.r_[np.zeros(n), -ub * np.ones(n), np.zeros(2 * n),
                 H @ w_bench - ind_dev, -(H @ w_bench + ind_dev), -tau]
    return Ain, bin_
obj = lambda x: 0.5 * x @ G @ x + dvec @ x
grad = lambda x: G @ x + dvec

# ---------------- 2. 自写凸 QP 原始-对偶内点法（16.9 节 + Mehrotra） ----------------
def qp_ipm(G, d, Aeq, beq, Ain, bin_, tol=1e-10, max_iter=100):
    n_, me, mi = len(d), len(beq), len(bin_)
    x = np.zeros(n_); y = np.maximum(Ain @ x - bin_, 1.0); lam = np.ones(mi); nu = np.zeros(me)
    for k in range(max_iter):
        rd = G @ x + d - Ain.T @ lam - Aeq.T @ nu
        re, ri = Aeq @ x - beq, Ain @ x - y - bin_
        mu = y @ lam / mi
        if max(np.abs(rd).max(), np.abs(re).max(), np.abs(ri).max(), mu) < tol:
            return x, lam, nu, k
        Dm = lam / y
        Kmat = np.block([[G + Ain.T @ (Dm[:, None] * Ain), -Aeq.T],
                         [Aeq, np.zeros((me, me))]])
        def solve(rc):
            rhs = np.r_[-rd - Ain.T @ ((rc + lam * ri) / y), -re]
            sol = np.linalg.solve(Kmat, rhs)
            dx, dnu = sol[:n_], sol[n_:]
            dy = Ain @ dx + ri
            dl = -(rc + lam * dy) / y
            return dx, dy, dl, dnu
        def maxstep(v, dv):
            neg = dv < 0
            return min(1.0, np.min(-v[neg] / dv[neg])) if neg.any() else 1.0
        dx, dy, dl, dnu = solve(y * lam)                         # 预测步 σ=0
        a = min(maxstep(y, dy), maxstep(lam, dl))                # 原始/对偶同一步长
        sigma = ((y + a * dy) @ (lam + a * dl) / mi / mu) ** 3
        dx, dy, dl, dnu = solve(y * lam + dy * dl - sigma * mu)  # 校正 + 中心化
        a = min(1.0, 0.99 * min(maxstep(y, dy), maxstep(lam, dl)))
        x, y, lam, nu = x + a * dx, y + a * dy, lam + a * dl, nu + a * dnu
    raise RuntimeError("IPM not converged")

tau = 0.30
Ain, bin_ = ineq(tau)
t0 = time.perf_counter(); x_ipm, lam, nu, it = qp_ipm(G, dvec, Aeq, beq, Ain, bin_); t_ipm = time.perf_counter() - t0

# ---------------- 3. scipy：SLSQP 与 trust-constr ----------------
cons_slsqp = [{"type": "eq", "fun": lambda x: Aeq @ x - beq, "jac": lambda x: Aeq},
              {"type": "ineq", "fun": lambda x: Ain @ x - bin_, "jac": lambda x: Ain}]
x_start = np.r_[w0, np.zeros(2 * n)]                     # 不交易：可行起点
t0 = time.perf_counter()
r1 = minimize(obj, x_start, jac=grad, constraints=cons_slsqp, method="SLSQP",
              options={"maxiter": 500, "ftol": 1e-12})
t_sl = time.perf_counter() - t0
t0 = time.perf_counter()
r2 = minimize(obj, x_start, jac=grad, hess=lambda x: G, method="trust-constr",
              constraints=[LinearConstraint(Aeq, beq, beq), LinearConstraint(Ain, bin_, np.inf)],
              options={"gtol": 1e-12, "xtol": 1e-14, "maxiter": 5000})
t_tc = time.perf_counter() - t0
print(f"{'方法':<14}{'目标值':>14}{'迭代':>6}{'用时s':>8}")
print(f"{'自写 IPM':<14}{obj(x_ipm):14.8f}{it:6d}{t_ipm:8.3f}")
print(f"{'SLSQP':<14}{r1.fun:14.8f}{r1.nit:6d}{t_sl:8.3f}")
print(f"{'trust-constr':<14}{r2.fun:14.8f}{r2.nit:6d}{t_tc:8.3f}")

w = x_ipm[:n]; turn = np.abs(w - w0).sum()
print(f"\n持仓数 {int((w > 1e-6).sum())}/{n}，顶格 {int((w > ub - 1e-6).sum())} 只，"
      f"换手 {turn:.4f}，预期 alpha {alpha @ w:.4%}，年化波动 {np.sqrt(w @ Sigma @ w):.2%}")
print("行业权重 - 基准:", (H @ w - H @ w_bench).round(4))
lam_turn = lam[-1]
print(f"换手约束乘子 λ_turn = {lam_turn:.5f}（每多 1 单位换手额度，目标改进量）")
x_b, *_ = qp_ipm(G, dvec, Aeq, beq, *ineq(tau + 1e-4))
print(f"有限差分验证: {(obj(x_ipm) - obj(x_b)) / 1e-4:.5f}")

print("trust-constr 状态:", r2.message, " 约束违反:", f"{r2.constr_violation:.1e}")

# ---------------- 4. 换手预算扫描 + 热启动 ----------------
from scipy.optimize import linprog
Ain, bin_ = ineq(10.0)
lp = linprog(np.r_[np.zeros(n), np.ones(2 * n)], A_ub=-Ain[:-1], b_ub=-bin_[:-1],
             A_eq=Aeq, b_eq=beq, bounds=[(None, None)] * N, method="highs")
print(f"\n满足行业/上限约束所需的最小换手 = {lp.fun:.4f}（低于它问题不可行）")
print("\n tau   效用(-obj)  换手乘子  SLSQP冷启动迭代  热启动迭代")
x_prev = x_start
for tau in [0.22, 0.30, 0.50, 0.80, 1.20]:
    Ain, bin_ = ineq(tau)
    cons = [{"type": "eq", "fun": lambda x: Aeq @ x - beq, "jac": lambda x: Aeq},
            {"type": "ineq", "fun": lambda x, A=Ain, b=bin_: A @ x - b, "jac": lambda x, A=Ain: A}]
    cold = minimize(obj, x_start, jac=grad, constraints=cons, method="SLSQP", options={"maxiter": 500, "ftol": 1e-12})
    hot = minimize(obj, x_prev, jac=grad, constraints=cons, method="SLSQP", options={"maxiter": 500, "ftol": 1e-12})
    xi, li, _, _ = qp_ipm(G, dvec, Aeq, beq, Ain, bin_)
    x_prev = hot.x
    print(f"{tau:5.2f} {-obj(xi):10.5f} {li[-1]:9.5f} {cold.nit:12d} {hot.nit:12d}")

输出：

方法                       目标值    迭代     用时s
自写 IPM            0.02172576    13   0.004
SLSQP             0.02172576    11   0.015
trust-constr      0.02172576    46   0.075

持仓数 34/40，顶格 1 只，换手 0.3000，预期 alpha 0.4697%，年化波动 11.36%
行业权重 - 基准: [-0.02 -0.02  0.02  0.02]
换手约束乘子 λ_turn = 0.09165（每多 1 单位换手额度，目标改进量）
有限差分验证: 0.09163
trust-constr 状态: `gtol` termination condition is satisfied.  约束违反: 6.9e-18

满足行业/上限约束所需的最小换手 = 0.1956（低于它问题不可行）

 tau   效用(-obj)  换手乘子  SLSQP冷启动迭代  热启动迭代
 0.22   -0.03060   0.13229            4            4
 0.30   -0.02173   0.09165           11           10
 0.50   -0.00620   0.06334            8            8
 0.80    0.00875   0.03886           23           21
 1.20    0.01845   0.00000           11            8

逐条解读：

三种方法目标值一致到第 8 位。自写内点法 13 次迭代，与第 14 章 LP 的表现相同；trust-constr 的障碍法需要 46 次；SLSQP 11 次。用时都在毫秒级（具体数字随机器而异），40 只股票的问题对任何方法都不难。规模到几千只、约束到几百条时，差别才会显现。
约束诊断。4 个行业的偏离全部顶在 ±2% 上，换手恰好用满 0.30——说明行业约束和换手预算都是有效约束，正在"限制"组合。只有 1 只股票顶格 8%，个股上限基本不是瓶颈。
换手的影子价格。乘子 \(\lambda_{\rm turn}=0.0917\)，有限差分 \(0.0916\)，两者一致：每多给 1 单位（100%）换手额度，目标（效用）改善约 0.092，即每多 1% 换手改善约 9bp 效用。这个数字可以和单边成本 20bp 比较：内点法已把 20bp 的成本计入目标，影子价格是扣除成本后的净边际收益。
换手–效用前沿。扫描 \(\tau\) 得到一条凹的前沿：影子价格随 \(\tau\) 增大单调下降（0.132 → 0.092 → 0.063 → 0.039 → 0），即边际收益递减；\(\tau=1.2\) 时换手约束不再有效，乘子为 0（互补松弛）。最小可行换手 0.1956 来自 LP，\(\tau\) 低于它时任何求解器都会报告不可行。
热启动。以前一个 \(\tau\) 的解作为初值，SLSQP 的迭代略有减少（23→21，11→8），但不显著。原因是 scipy 的 SLSQP 并不保存工作集和分解，只是从更好的点出发；真正保留工作集的有效集 QP 求解器从热启动中获益要大得多（16.5.7 节）。

金融直觉：把这个模型的 KKT 条件逐股写出来，可以得到一个实务上非常有用的结论：不交易区间。记预算乘子 \(\nu\)、换手拆分等式 \(w_i-b_i+s_i=w_i^0\) 的乘子 \(\eta_i\)、换手预算的乘子 \(\lambda_T\)（按原书约定写成 \(\tau-\mathbf 1^T(b+s)\ge0\)），并定义股票 \(i\) 的边际净吸引力 \(m_i=\alpha_i-2\kappa(\Sigma w)_i+\nu+(\text{行业约束乘子带来的项})\)，即"多持有一点这只股票带来的 alpha，减去它增加的风险，再按预算和行业约束的影子价格调整"。对权重严格介于 0 和上限之间的股票（顶格或清仓的股票还要再加上界约束乘子，逻辑相同），由 \(w_i\) 的平稳性得 \(m_i=-\eta_i\)，再由 \(b_i\)、\(s_i\) 的平稳性（\(c+\eta_i+\lambda_T=\beta_i\)、\(c-\eta_i+\lambda_T=\zeta_i\)，\(\beta_i,\zeta_i\ge0\) 是 \(b_i\ge0\)、\(s_i\ge0\) 的乘子）配合互补松弛得：买入的股票（\(b_i>0\)）：\(m_i=c+\lambda_T\)；卖出的股票（\(s_i>0\)）：\(m_i=-(c+\lambda_T)\)；不交易的股票：\(-(c+\lambda_T)\le m_i\le c+\lambda_T\)。读法：只有边际净吸引力超过"单边成本 + 换手预算的影子价格"，才值得交易；在这个区间之内，哪怕模型认为它略被高估或低估，也按兵不动。所以换手预算在效果上等于把交易成本从 \(c\) 抬高到 \(c+\lambda_T\)。本例中 \(c=0.002\)、\(\lambda_T\approx0.092\)，换手预算隐含的"成本"约为显性费率的 46 倍（都以目标函数的效用单位计），说明在 \(\tau=0.30\) 时真正卡住组合的是换手额度，而不是佣金和冲击成本。这也解释了上面第 4 条：放宽换手预算，\(\lambda_T\) 下降，不交易区间收窄，更多股票被调到目标附近。这个结构和第 12 章 12.10.3 节的等边际原则表一脉相承，只是"边际相等"被放宽成了"边际落在一个带宽内"。

16.11.3 代码二：梯度投影法求解多空界约束组合

市场中性的多空组合常把净敞口用期货对冲，个股只受多空限额约束：\(\min\kappa w^T\Sigma w-\alpha^Tw\) s.t. \(-2\%\le w_i\le2\%\)。这是标准的界约束 QP (16.43)。

import numpy as np
from scipy.optimize import minimize

def cauchy_point(G, d, x, l, u):
    """沿投影梯度路径 x(t)=P(x - t g) 找第一个局部极小点 (16.44)-(16.48)"""
    g = G @ x + d
    with np.errstate(divide="ignore", invalid="ignore"):
        tbar = np.where(g < 0, (x - u) / g, np.where(g > 0, (x - l) / g, np.inf))
    tbar = np.where(np.isnan(tbar), np.inf, tbar)
    bps = np.unique(tbar[(tbar > 0) & np.isfinite(tbar)])
    t_prev, xc = 0.0, x.copy()
    for t_next in np.r_[bps, np.inf]:
        p = np.where(tbar > t_prev, -g, 0.0)            # 本段方向 (16.47)
        f1 = d @ p + xc @ G @ p                          # f'
        f2 = p @ G @ p                                   # f''
        if f1 >= 0:                                      # 本段起点即极小
            return xc
        dt = -f1 / f2 if f2 > 0 else np.inf
        if dt < t_next - t_prev:
            return xc + dt * p
        if not np.isfinite(t_next):
            raise ValueError("unbounded")
        xc = np.clip(xc + (t_next - t_prev) * p, l, u)
        t_prev = t_next
    return xc

def gradient_projection(G, d, l, u, x0, tol=1e-10, max_iter=200):
    """算法 16.2：Cauchy 点 + 自由变量上的 CG 子空间极小化（碰界即停）"""
    x = np.clip(x0, l, u)
    for k in range(max_iter):
        g = G @ x + d
        pg = np.where((x <= l) & (g > 0), 0, np.where((x >= u) & (g < 0), 0, g))
        if np.linalg.norm(pg, np.inf) < tol:
            return x, k
        xc = cauchy_point(G, d, x, l, u)
        free = (xc > l + 1e-12) & (xc < u - 1e-12)
        # 在自由变量上从 xc 出发做 CG，求 (16.49) 的近似解
        z = xc.copy(); r = (G @ z + d)[free]; p = -r
        for _ in range(free.sum()):
            if np.linalg.norm(r) < 1e-14: break
            Gp = G[np.ix_(free, free)] @ p
            curv = p @ Gp
            a = r @ r / curv
            zf = z[free] + a * p
            if np.any(zf < l[free]) or np.any(zf > u[free]):   # 碰界：截断到可行并停止
                lo = np.where(p < 0, (l[free] - z[free]) / p, np.inf)
                hi = np.where(p > 0, (u[free] - z[free]) / p, np.inf)
                a = min(a, np.min(np.minimum(lo, hi)))
                z[free] = z[free] + a * p
                break
            z[free] = zf
            r_new = r + a * Gp
            p = -r_new + (r_new @ r_new) / (r @ r) * p
            r = r_new
        x = np.clip(z, l, u)
    return x, max_iter

rng = np.random.default_rng(11)
n = 300
B = rng.standard_normal((n, 5)) * 0.2
Sigma = B @ B.T + np.diag(rng.uniform(0.02, 0.08, n))
alpha = 0.03 * rng.standard_normal(n)
kappa = 5.0
G, d = 2 * kappa * Sigma, -alpha                     # min κ w'Σw - α'w
l, u = -0.02 * np.ones(n), 0.02 * np.ones(n)         # 多空个股限额 ±2%
x, it = gradient_projection(G, d, l, u, np.zeros(n))
print(f"梯度投影: 外迭代 {it} 次，目标 {0.5*x@G@x+d@x:.8f}，"
      f"顶格 {int((np.abs(x) > 0.02 - 1e-9).sum())}/{n} 只")
r = minimize(lambda w: 0.5*w@G@w + d@w, np.zeros(n), jac=lambda w: G@w + d,
             method="L-BFGS-B", bounds=list(zip(l, u)), options={"ftol": 1e-15, "gtol": 1e-12})
print(f"L-BFGS-B: 迭代 {r.nit} 次，目标 {r.fun:.8f}，与梯度投影解差 {np.abs(r.x-x).max():.1e}")

输出：

梯度投影: 外迭代 41 次，目标 -0.11572584，顶格 210/300 只
L-BFGS-B: 迭代 114 次，目标 -0.11572584，与梯度投影解差 1.1e-07

300 只股票中有 210 只顶在 ±2% 的限额上。经典有效集法从 \(w=0\) 出发每步只能加一个约束，至少需要 210 次迭代；梯度投影法只用了 41 次外迭代，因为每次 Cauchy 点计算都能一次把许多变量推到界上。L-BFGS-B（原书第 9 章）也是"投影 + 子空间"的思路，只是子空间里用有限内存拟牛顿近似代替精确 Hessian，结果一致。

16.11.4 实务清单

先检查凸性：协方差矩阵是否半正定（最小特征值），否则先修正（16.6 节）。
先检查可行性：约束多时用 LP 求最小换手、最小跟踪误差等，确认问题可行再优化。
缩放：权重量级 \(10^{-2}\)、方差量级 \(10^{-2}\)–\(10^{-4}\)、成本量级 \(10^{-3}\)，乘子的大小依赖缩放（第 16a 章练习 5）。把收益和风险统一成年化或日度、把权重用百分比，能显著改善求解器的稳定性。
选求解器：中小规模、日频再平衡、可热启动 → 有效集法；大规模、多约束、带锥约束 → 内点法；只有界约束 → 梯度投影 / L-BFGS-B；超大规模可接受中等精度 → 一阶方法（OSQP 的 ADMM，即第 17 章增广拉格朗日的分块版本）。
读乘子：每条约束的乘子就是它的影子价格，是和投资经理、风控讨论约束松紧时最有用的数字。

本章（下）小结

梯度投影法沿投影最速下降的分段路径找 Cauchy 点，再在自由变量上做子空间极小化，一次迭代就能让大量约束进入或离开有效集，最适合界约束 QP（也适合严格凸 QP 的对偶）。凸 QP 的原始–对偶内点法是 LP 内点法的直接推广：KKT 条件加松弛 \(y\)，朝中心路径 \(y_i\lambda_i=\sigma\mu\) 走牛顿步，每步解增广系统或正规方程 \((G+A^TY^{-1}\Lambda A)\Delta x=\cdots\)；与 LP 不同，原始与对偶必须用同一步长。有效集法步多而便宜、易热启动，内点法步少而贵、结构固定、适合大规模。严格凸 QP 的对偶是 \(\lambda\ge0\) 的界约束 QP。组合优化实例显示：含换手拆分、行业偏离、个股上限的均值–方差问题可以被 SLSQP、trust-constr 和自写内点法一致地求解，换手约束的乘子精确等于效用对换手预算的边际改进，随预算增加单调递减直至为零。

概念	公式 / 要点
投影	\(P(x,l,u)_i=\mathrm{mid}(l_i,x_i,u_i)\)
投影路径	\(x(t)=P(x-tg,l,u)\)，断点 \(\bar t_i\) (16.46)
Cauchy 点	路径上第一个局部极小：\(f'\ge0\) 停，或段内 \(\Delta t^*=-f'/f''\)
子空间极小化	固定 \(\mathcal A(x^c)\)，自由变量上 CG，碰界即停；只需 \(q(x^+)\le q(x^c)\)
QP 的 KKT（松弛形式）	\(Gx-A^T\lambda+d=0\)，\(Ax-y-b=0\)，\(y_i\lambda_i=0\)，\((y,\lambda)\ge0\)
内点步	(16.53)；正规方程 \((G+A^TY^{-1}\Lambda A)\Delta x=\cdots\)
QP 内点的步长	原始、对偶同一步长（\(G\) 耦合）
QP 对偶	\(\max_{\lambda\ge0}-\frac12\lambda^TAG^{-1}A^T\lambda+\lambda^T(b+AG^{-1}d)-\frac12d^TG^{-1}d\)
换手拆分	\(w-w^0=b-s\)，\(b,s\ge0\)，成本 \(>0\) 时 \(b_is_i=0\)

练习

基础

证明：界约束问题的 KKT 条件等价于投影梯度为零：\(P(x-g,l,u)=x\)。
对二维问题 \(\min\frac12(x_1^2+x_2^2)-3x_1-x_2\) s.t. \(0\le x\le1\)，从 \(x^0=(0,0)\) 出发手算 Cauchy 点，并判断它是否已是最优解。提示：\(g=(-3,-1)\)，\(\bar t=(1/3,1)\)；第一段 \(f'=-10<0\)，\(f''=10\)，\(\Delta t^*=1>1/3\)，进入第二段。
写出 16.9 节对等式 + 不等式混合约束凸 QP 的 KKT 条件，并推导类似 (16.53) 的原始–对偶步方程（与 16.11.2 中 qp_ipm 的实现对照）。
由 (16.56) 消去 \(x\) 推导 (16.57)，并说明为什么需要 \(G\) 正定。
在 16.11.2 中把单边费率改为 0，观察最优解是否仍满足 \(b_is_i=0\)；如果不满足，换手约束的含义会发生什么变化？

进阶

对"长仓 + 预算"约束集 \(\{w\ge0,\ \mathbf 1^Tw=1\}\)（概率单纯形），推导欧氏投影的排序算法：存在阈值 \(\theta\) 使 \(P(v)_i=\max(v_i-\theta,0)\)。用它实现投影梯度法求解长仓最小方差组合，并与 SLSQP 比较。
在 16.11.2 的模型中加入跟踪误差约束 \((w-w_b)^T\Sigma(w-w_b)\le\mathrm{TE}^2\)。它不再是 QP，而是二次约束问题（QCQP/SOCP）。用 trust-constr 或 SLSQP 求解，画出跟踪误差上限与预期 alpha 的前沿，并读出该约束的乘子。
用 16.10 节的对偶方法求解无界约束、只有 \(K\) 条因子中性约束 \(B^Tw\ge-\epsilon\)、\(-B^Tw\ge-\epsilon\) 的均值–方差问题。比较原始问题（\(n=2000\)）与对偶问题（\(2K=20\) 维）的求解时间。
修改 qp_ipm，把每步的稠密求解换成"因子模型 + Woodbury"的结构化求解：利用 \(G=2\kappa(BFB^T+D)\)，推导正规方程系数矩阵的低秩加对角结构。
把 16.11.3 的问题改为"长仓 + 预算 + 个股上限"，比较三种做法：(a) 用罚项 \(\frac{\rho}2(\mathbf 1^Tw-1)^2\) 处理预算再用 L-BFGS-B；(b) 用增广拉格朗日（第 17 章）；(c) 用 SLSQP。观察 (a) 中 \(\rho\) 增大时约束违反和迭代次数的变化。

原书推荐习题：16.17（求 (16.49) 的零空间基）、16.18（混合约束凸 QP 的原始–对偶内点步，组合优化求解器的核心）、16.12（编程实现有效集法求解给定 QP）、16.14（乘子与缩放）。

原书对照

本章内容	原书位置	PDF 页码
16.8 梯度投影法、Cauchy 点、子空间极小化、算法 16.2	§16.6 The Gradient-Projection Method，式 (16.43)–(16.49)	PDF p.493–498
16.9 凸 QP 内点法、增广系统与正规方程	§16.7 Interior-Point Methods，式 (16.50)–(16.54)	PDF p.498–501
16.10 QP 对偶	§16.8 Duality，式 (16.55)–(16.57)	PDF p.501–502
注释与参考（NP-hard、Markowitz、带状结构等）	Notes and References	PDF p.502–503
习题 16.1–16.18	Exercises	PDF p.503–505

说明：原书 (16.53) 在抽取文本中矩阵列序错乱，本章按含义还原；习题 16.16 原文表述（"\(\bar Z=[Z,z]\) 时 \(\bar Z^TW\bar Z\) 正定"）一般不成立，应理解为删去一列（即加约束）时保持正定。原书正文页码 = PDF 页码减 19（已用 PDF 页眉核对，第 12–16 章均如此）。

第 16b 章 二次规划：梯度投影法、内点法、对偶与组合优化实战