第 18 章 序列二次规划

学习目标

读完本章，你应当能够：

1. 从"对 KKT 方程组用牛顿法"推导局部 SQP，说明它与"在线性化约束下极小化 Lagrange 函数的二次模型"等价，并知道它局部二次收敛。
2. 把 SQP 推广到不等式约束（IQP 与 EQP 两种思路），理解接近解时有效集会稳定下来（定理 18.1）。
3. 掌握步计算的三种方式（增广系统、值空间、零空间）、最小二乘乘子和约化 Hessian 方法。
4. 知道 Hessian 的几种选择：精确 Lagrange Hessian、阻尼 BFGS、SR1、增广拉格朗日 Hessian、约化 Hessian 拟牛顿；理解阻尼 BFGS 如何保证正定。
5. 用 \(\ell_1\) 或 Fletcher 价值函数做线搜索，会按 \(1/\mu>\|\lambda\|_\infty\) 更新罚参数；了解信赖域 SQP 的三种处理线性化约束不相容的办法（约束平移、双椭球、S\(\ell_1\)QP）。
6. 理解 Maratos 效应及其三种对策（Fletcher 价值函数、二阶校正、watchdog），并能读懂 scipy.optimize.minimize(method="SLSQP") 的行为与报错。

读前导读

这一章在解决什么问题。 均值–方差优化是二次规划（QP），你在 CFA 里见过，求解器也很成熟。可一旦加入"组合波动 ≤ 12%""单资产风险贡献 ≤ 10%""冲击成本按 \(|\Delta w|^{3/2}\) 增长"这类约束，问题就不再是 QP。SQP 的办法是：在当前点把原问题近似成一个 QP（目标取二次近似、约束取一次近似），解这个 QP 得到一步，走过去再近似，如此反复。scipy 里最常用的 SLSQP 就是这个算法，所以本章也是在解释你日常调用的工具内部发生了什么、它的报错是什么意思。

和你熟悉的概念对照：久期是价格对收益率的一阶近似，加上凸性就是二阶近似；SQP 对目标做的正是"久期 + 凸性"式的二阶近似，对约束只做"久期"式的一阶近似。本章的很多麻烦（线性化约束不相容、Maratos 效应）都源于约束只做了一阶近似，而真实约束是弯的。

需要先想起来的数学。

• 牛顿法解方程组。 解 \(F(z)=0\) 时，在当前点用一阶泰勒展开 \(F(z+\Delta)\approx F(z)+J(z)\Delta\)，令其为零得 \(J\Delta=-F\)，\(J\) 是 \(F\) 的雅可比矩阵（各分量对各变量的偏导排成的矩阵）。一维例子：解 \(z^2-2=0\)，从 \(z=1.5\) 出发，\(\Delta=-(1.5^2-2)/(2\times1.5)\approx-0.083\)，一步到 \(1.4167\)，已经接近 \(\sqrt2=1.4142\)。离解足够近时误差每步平方（"二次收敛"，有效数字翻倍）。见 第 00 册第 02 章 导数与泰勒展开 和本册第 11 章。
• KKT 条件与拉格朗日函数。 本书记号 \(\mathcal L=f-\lambda^Tc\)，最优点满足 \(\nabla f=A^T\lambda\)（目标梯度是约束梯度的线性组合）和 \(c=0\)。乘子 \(\lambda\) 是影子价格。\(W=\nabla^2_{xx}\mathcal L\) 是拉格朗日函数对 \(x\) 的 Hessian。见 第 00 册第 05 章 多元微积分与优化 和本册第 12 章。
• 零空间与值空间。 约束雅可比 \(A\)（\(m\times n\)）的零空间 \(\{d:Ad=0\}\) 是"沿约束面移动"的方向，\(Z\) 的列是它的一组基；值空间方向（\(A^T\) 的列张成）是"垂直于约束面"的方向，\(Y\) 常取 \(A^T\)。任何步都能拆成 \(p=Yp_Y+Zp_Z\)。投影矩阵 \(P=I-A^T(AA^T)^{-1}A\) 把向量投到零空间。见 第 00 册第 06 章 线性代数速成 和第 16a 章。
• 向量范数。 \(\|c\|_1=\sum|c_i|\)，\(\|c\|_2=\sqrt{\sum c_i^2}\)，\(\|\lambda\|_\infty=\max|\lambda_i|\)。Hölder 不等式 \(|c^T\lambda|\le\|c\|_1\|\lambda\|_\infty\)：内积不超过"总量 × 最大单价"。例如 \(c=(1,-2)\)、\(\lambda=(3,1)\)，\(|c^T\lambda|=1\le3\times3=9\)。
• BFGS 与曲率条件。 拟牛顿法用梯度差 \(y\) 和步 \(s\) 更新 Hessian 近似 \(B\)；只有 \(s^Ty>0\) 时 BFGS 才能保持 \(B\) 正定。见本册第 08 章。

怎么读这一章。 核心必读是 18.1（SQP = 对 KKT 用牛顿法，这是全章的地基）、18.4.1（阻尼 BFGS，SLSQP 用的就是它）、18.5.1（\(\ell_1\) 价值函数和"罚权重要超过最大乘子"）、18.10（Maratos 效应），再加 18.11 的两段代码和解读。第一次可以只看结论：18.3.1 的三种解法细节、18.4.3 与 18.7 的约化 Hessian 方法、18.5.2 Fletcher 函数的方向导数 (18.38)、18.8 信赖域 SQP 的三种重构、18.9 的收敛速度定理。建议顺序：18.1 → 18.2 → 18.4.1 → 18.5.1 → 18.6 → 18.10 → 18.11，其余按需回看。

18.1 局部 SQP 方法

先说结论：序列二次规划（sequential quadratic programming, SQP）在每个迭代点用一个二次规划子问题近似原问题，以子问题的解作为步。它是非线性约束优化最有效的方法之一，可以放在线搜索或信赖域框架中，适用于小规模和大规模问题。与上一章的 SLC 方法（多数约束线性时有效）不同，SQP 在约束显著非线性时更有优势。

18.1.1 从牛顿法出发

先看等式约束问题

\(f:\mathbb R^n\to\mathbb R\)，\(c:\mathbb R^n\to\mathbb R^m\) 光滑。Lagrange 函数 \(\mathcal L(x,\lambda)=f(x)-\lambda^Tc(x)\)，约束 Jacobian \(A(x)^T=[\nabla c_1(x),\dots,\nabla c_m(x)]\)。KKT 条件是 \(n+m\) 个方程、\(n+m\) 个未知数：

对它用牛顿法。\(F\) 的 Jacobian 为 \(\begin{bmatrix}W(x,\lambda)&-A(x)^T\\A(x)&0\end{bmatrix}\)，\(W=\nabla^2_{xx}\mathcal L\)，牛顿步 \((p_k,p_\lambda)\) 满足

推导拆解：(18.7) 就是"雅可比 × 步 = −残差"。把未知数看成 \(z=(x,\lambda)\)，对 \(F\) 的两块分别求导：

• 第一块 \(\nabla f(x)-A(x)^T\lambda=\nabla f-\sum_i\lambda_i\nabla c_i\)：对 \(x\) 求导得 \(\nabla^2f-\sum_i\lambda_i\nabla^2c_i=W\)；对 \(\lambda\) 求导得 \(-A^T\)（它对 \(\lambda\) 是线性的）。
• 第二块 \(c(x)\)：对 \(x\) 求导得 \(A\)；与 \(\lambda\) 无关，导数为 \(0\)。

于是雅可比为 \(\begin{bmatrix}W&-A^T\\A&0\end{bmatrix}\)，右端是 \(-F=\begin{bmatrix}-\nabla f+A^T\lambda\\-c\end{bmatrix}\)。注意 \(W\) 里含有约束的二阶导 \(\nabla^2c_i\)，权重是乘子 \(\lambda_i\)——这是后面"约束曲率通过乘子进入模型"的来源。

这称为 Newton–Lagrange 方法。KKT 矩阵非奇异的条件是：

假设 18.1：(a) \(A_k\) 行满秩（LICQ）；(b) \(W_k\) 在约束切空间上正定：\(d^TW_kd>0\)，\(\forall d\ne0,\ A_kd=0\)。

（与第 16a 章引理 16.1 完全相同。）在满足二阶充分条件的解附近这两条成立，牛顿迭代二次收敛。

18.1.2 SQP 的观点

在 \((x_k,\lambda_k)\) 处定义 QP 子问题

在假设 18.1 下它有唯一解 \((p_k,\mu_k)\)（\(\mu_k\) 是子问题的乘子），满足 \(W_kp_k+\nabla f_k-A_k^T\mu_k=0\)，\(A_kp_k+c_k=0\)。把 (18.7) 第一行两边减去 \(A_k^T\lambda_k\)，并令 \(\lambda_{k+1}=\lambda_k+p_\lambda\)：

这正是 (18.8) 的 KKT 系统。所以 SQP 与对 KKT 条件的牛顿法等价：新的 \(x\) 是 QP 的解，新的 \(\lambda\) 是 QP 的乘子。牛顿观点便于分析，SQP 观点便于设计实用算法并推广到不等式。

另一种理解：(18.8) 的线性项可以换成 \(\nabla_x\mathcal L(x_k,\lambda_k)^Tp\)（在约束下二者只差常数），此时目标就是 Lagrange 函数的二次近似——SQP 是在线性化约束下极小化 Lagrange 函数的二次模型。注意二次项用的是 Lagrange 函数的 Hessian \(W_k\)，而不是 \(\nabla^2f\)：约束的曲率通过乘子进入了模型。若只用 \(\nabla^2f\)，算法会在约束弯曲时失去快速收敛。

推导拆解：从 (18.7) 到 (18.10) 只有一步代数。(18.7) 第一行是 \(W_kp_k-A_k^Tp_\lambda=-\nabla f_k+A_k^T\lambda_k\)，把右端的 \(A_k^T\lambda_k\) 移到左边：\(W_kp_k-A_k^T(\lambda_k+p_\lambda)=-\nabla f_k\)，再记 \(\lambda_{k+1}=\lambda_k+p_\lambda\) 即得。另一方面，QP (18.8) 的拉格朗日函数是 \(\frac12p^TW_kp+\nabla f_k^Tp-\mu^T(A_kp+c_k)\)，对 \(p\) 求导令其为零得 \(W_kp+\nabla f_k-A_k^T\mu=0\)，再加上约束 \(A_kp+c_k=0\)，正好是 (18.10)。两边方程一模一样，所以解也一样。

金融直觉：用例 18.1（单位圆约束）看 \(W_k\) 与 \(\nabla^2f\) 的差别。\(f=2(x_1^2+x_2^2-1)-x_1\)，\(\nabla^2f=4I\)；约束 \(c=x_1^2+x_2^2-1\)，\(\nabla^2c=2I\)；\(\lambda^*=1.5\)，所以 \(W=4I-1.5\times2I=I\)。只看 \(\nabla^2f\) 会把曲率高估 4 倍，牛顿步会短很多。这好比给一个带利率约束的组合算对冲比率：只看资产本身的凸性、忽略约束（如久期匹配要求）本身的凸性，二阶调整就会算错。乘子 \(\lambda\) 相当于约束的"头寸规模"，约束越值钱，它的曲率在模型里的权重越大。

算法 18.1（局部 SQP）

选初始 (x0, λ0)；
for k = 0,1,2,...
    计算 fk, ∇fk, Wk = W(xk, λk), ck, Ak；
    解 (18.8) 得 pk 与乘子 μk；
    x_{k+1} = xk + pk；λ_{k+1} = μk；
    若满足收敛测试：STOP；
end

若解处假设 18.1 成立、二阶导 Lipschitz 连续、初始点足够近，则 \((x_k,\lambda_k)\) 二次收敛（定理 18.4，由牛顿法理论直接得到）。

18.1.3 不等式约束

一般问题

的子问题线性化所有约束：

用第 16a、16b 章的 QP 算法求解，\(p_k\) 与 \(\lambda_{k+1}\) 取其解与乘子。

定理 18.1：设 \(x^*\) 是 (18.11) 的解，有效约束 Jacobian \(A_*\) 行满秩，\(W_*\) 在 \(A_*\) 的零空间上正定，严格互补成立。则 \((x_k,\lambda_k)\) 足够接近 \((x^*,\lambda^*)\) 时，子问题 (18.12) 有一个局部解，其有效集与 \(\mathcal A(x^*)\) 相同。

即接近解时有效集固定，SQP 的行为与等式约束情形相同，前面的收敛理论直接适用。

两种实现思路：

• IQP（不等式约束 QP）：每步解完整的 (18.12)，以其解的有效集作为最优有效集的猜测。实践中很成功；缺点是大问题中解一般 QP 代价高，但接近解时用上一步信息热启动会很便宜。
• EQP（等式约束 QP）：每步先选一个工作集，只解工作集约束作等式的子问题 (18.8)，其余约束暂时忽略；工作集用乘子符号或辅助子问题更新。例如 SLP-EQP：先去掉二次项、加信赖域，解一个 LP 确定工作集，再固定工作集、加回二次项解等式 QP。

18.2 实用 SQP 方法概览

局部 SQP 只在解附近有效。实用 SQP 还要能从远处起点出发、处理非凸问题。类比无约束优化：牛顿模型在解附近 Hessian 正定时好用，远离解时可能非凸；信赖域法限制步长，线搜索法把 Hessian 修正为正定（或用拟牛顿近似）。

• 线搜索 SQP：\(W_k\) 在约束切空间上不正定时，用正定近似 \(B_k\) 代替，或在分解时直接修改，或改用具有凸性的增广拉格朗日 Hessian。
• 信赖域 SQP：子问题加信赖域约束，可以直接用不定的 \(W_k\)；但信赖域可能使线性化约束不相容，需要放松约束，增加复杂性。两类方法各有取舍。
• 价值函数：用第 15 章的价值函数判断步的好坏；参数（罚参数 \(\mu\)）的更新规则对实际表现影响很大。

18.3 步的计算

18.3.1 等式约束

解 KKT 系统 (18.10)（\(W_k\) 可以是精确 Hessian 或拟牛顿近似 \(B_k\)）。三种方式与第 16a 章 16.3 节一一对应：

• 增广系统法：对 KKT 矩阵做对称不定分解 \(LDL^T\)（稠密用 Bunch–Kaufman，稀疏用 Duff–Reid），直接得到 \(p_k\) 和 \(\lambda_{k+1}\)。
• 值空间法：\(W_k\) 正定时

对显式维护逆 Hessian 近似 \(H_k\approx W_k^{-1}\) 的拟牛顿法特别方便。

• 零空间法（许多 SQP 软件的核心）：\(p_k=Y_kp_Y+Z_kp_Z\)，

乘子 \((A_kY_k)^T\lambda_{k+1}=Y_k^T(\nabla f_k+W_kp_k)\)（18.15）。只要求约化 Hessian \(Z_k^TW_kZ_k\) 正定。

金融直觉：以预算约束 \(\mathbf 1^Tw=1\) 为例，\(A=\mathbf 1^T\)。零空间方向 \(Z\) 是"加总为零"的调仓——买一只、卖另一只，总投入不变，相当于自融资交易；值空间方向 \(Y=A^T=\mathbf 1\) 是"所有仓位同比例加减"，只改变总投入。零空间法的两步由此很好理解：先用 \(p_Y\) 把总投入拉回到 1（修正约束违反，\((AY)p_Y=-c\)）；再在所有自融资交易中挑出最能改善目标的那一个（\(p_Z\)，只看目标在约束面内的曲率 \(Z^TWZ\)）。\(Z^TWZ\) 正定，意思是"任何自融资调仓都会让二次模型变差"在最优点成立，这比要求 \(W\) 本身正定宽松得多。

白话解释：(18.16) 的"最小二乘乘子"在问：用约束梯度的线性组合 \(A^T\lambda\) 去拟合目标梯度 \(\nabla f\)，最优系数是多少？这和 OLS 回归 \(\hat\beta=(X^TX)^{-1}X^Ty\) 是同一个公式，\(X=A^T\)、\(y=\nabla f\)。在最优点拟合残差为零（KKT 条件），离最优点越远残差越大。

最小二乘乘子。接近解时 \(p_k\to0\) 而 \(\nabla f_k\) 一般不趋于零，所以可以去掉 (18.15) 右端的 \(W_kp_k\)，取 \(Y_k=A_k^T\)，得

即 \(\min_\lambda\|\nabla f_k-A_k^T\lambda\|_2\) 的解。远离解时它也有意义（尽量满足一阶条件）。实践中在 \(x_{k+1}\) 处计算梯度后再用 (18.16) 求 \(\lambda_{k+1}\)，这样 SQP 就从 \((x,\lambda)\) 迭代变成纯原始的 \(x\) 迭代。

约化 Hessian 方法。再去掉交叉项 \(Z_k^TW_kY_kp_Y\)：

只需存储和近似 \((n-m)\times(n-m)\) 的约化 Hessian。合理性在于：法向分量 \(p_Y\) 通常比切向分量 \(p_Z\) 收敛得快（18.7 节）。

18.3.2 不等式约束与热启动

用有效集 QP 方法（算法 16.1；\(W_k\) 非凸时用不定 QP 变体）解 (18.12)。热启动对线搜索 SQP 的效率至关重要：用上一个子问题的解作初始点，用上一次 SQP 迭代的最终有效集作初始工作集，只有线性约束时还可以复用矩阵分解。

18.3.3 线性化约束不相容

线性化可能使子问题不可行。例：\(n=1\)，约束 \(x\le1\)、\(x^2\ge0\)，在 \(x_k=3\) 处线性化得 \(3+p\le1\) 与 \(9+6p\ge0\)，即 \(p\le-2\) 且 \(p\ge-1.5\)，矛盾——尽管原问题本身是可行的。

SNOPT 的弹性模式：若 (18.12) 不可行，改解

\(\gamma\) 为罚参数（\(\ell_1\) 精确罚的思想）。原问题可行且 \(\gamma\) 足够大时两者解相同；原问题不可行时，通常给出一个"违反最少"的点。SNOPT 取 \(\gamma=100\|\nabla f(x_s)\|\)，\(x_s\) 是首次检测到不相容的迭代点。

量化提示：scipy 的 SLSQP 遇到这种情况会报 "Inequality constraints incompatible"（status 4）或 "Positive directional derivative for linesearch"（status 8）。它没有弹性模式，所以要么换一个更可行的起点，要么把可能冲突的约束软化。

内点法解子问题：对一般问题，内点法只在早期迭代（有效集变化大、热启动收益小）与有效集法竞争，因为它不能像有效集法那样利用先验信息。

18.4 二次模型的 Hessian

取 \(W_k=\nabla^2_{xx}\mathcal L(x_k,\lambda_k)\) 能获得二次收敛，远离解时也常进展迅速；但需要二阶导（可能难算），而且在约束零空间上不一定正定。替代选择如下。

18.4.1 完整拟牛顿近似与阻尼 BFGS

对 \(\nabla^2_{xx}\mathcal L\) 维护 BFGS 近似 \(B_k\)，

（注意两个梯度用同一个 \(\lambda_{k+1}\)，相当于把乘子固定的 Lagrange 函数当作目标。）若 \(\nabla^2_{xx}\mathcal L\) 在相关区域正定，效果如同无约束 BFGS。但 Lagrange Hessian 常常不正定（它只需在切空间上正定），此时曲率条件 \(s_k^Ty_k>0\) 即使在解附近也可能不成立，标准 BFGS 无法保持正定。

• 跳过更新：若 \(s_k^Ty_k<\theta s_k^TB_ks_k\)（如 \(\theta=10^{-2}\)）就跳过。在许多问题上可行，但不够通用。
• 阻尼 BFGS（Powell）：用 \(y_k\) 与 \(B_ks_k\) 的插值代替 \(y_k\)：

推导拆解：为什么 \(\theta_k\ne1\) 时 \(s_k^Tr_k=0.2\,s_k^TB_ks_k\)？记 \(a=s_k^TB_ks_k>0\)（\(B_k\) 正定），\(b=s_k^Ty_k\)，此时 \(b<0.2a\)。

1. 内积是线性的：\(s_k^Tr_k=\theta_kb+(1-\theta_k)a=a-\theta_k(a-b)\)。
2. 代入 \(\theta_k=0.8a/(a-b)\)（因 \(b<0.2a<a\)，分母为正）：\(s_k^Tr_k=a-0.8a=0.2a\)。

所以 \(r_k\) 是在"真实曲率 \(y_k\)"和"旧模型曲率 \(B_ks_k\)"之间插值，权重刚好让曲率条件以 \(0.2a\) 的余量成立。数值例子：\(a=1\)、\(b=-0.5\)（负曲率），\(\theta=0.8/1.5\approx0.533\)，\(s^Tr=0.533\times(-0.5)+0.467\times1=0.2\)。可以把它理解成对一个噪声很大的波动率估计做收缩：新信息（\(y_k\)）不可信时，向先验（\(B_ks_k\)）收缩到刚好安全的程度。

\(\theta_k\ne1\) 时恰有 \(s_k^Tr_k=0.2\,s_k^TB_ks_k>0\)（18.24），所以 \(B_{k+1}\) 一定正定。\(\theta_k=0\) 时 \(B_{k+1}=B_k\)，\(\theta_k=1\) 是原 BFGS。这就是 SLSQP 等许多 SQP 程序采用的策略。缺点是没有解决 Lagrange Hessian 本身不正定的根本问题，困难问题上仍可能表现很差。此时 SR1 更合适（它允许不定近似），是信赖域 SQP 的好选择。

18.4.2 增广拉格朗日 Hessian

由第 17 章定理 17.5，

在 \(\mu\) 小于某阈值时正定：罚项在 \(A^T\) 的列空间方向上加了正曲率，零空间上的曲率不变（而零空间上本来就正定）。所以可以用 \(\nabla^2\mathcal L_A\) 或其拟牛顿近似作 \(W_k\)。困难是阈值依赖未知量：\(\mu\) 太小时罚项主导、表现差；太大时可能不正定。

18.4.3 约化 Hessian 拟牛顿近似

只近似 \((n-m)\) 维的 \(Z_k^T\nabla^2_{xx}\mathcal L Z_k\)——在标准假设下它在解附近正定，所以 BFGS 很自然。迭代为

\(M_k\) 是约化 Hessian 的近似，用割线方程 \(M_{k+1}s_k=y_k\) 和 BFGS 更新：

Coleman–Conn 变体只沿约束切空间收集曲率：

需要在中间点多求一次梯度，但可以证明解附近 \(y_k^Ts_k>0\)，可以安全地做 BFGS。

18.5 价值函数与下降性

价值函数 \(\phi\) 控制步长（线搜索）或决定是否接受步（信赖域），扮演无约束优化中目标函数的角色。下面限于等式约束。

18.5.1 \(\ell_1\) 价值函数

它在 \(c_i(x)=0\) 处不可微，但总有方向导数（附录 A 给出 \(\|x\|_1\) 的方向导数公式）。

引理 18.2：若 \(p_k,\lambda_{k+1}\) 由 (18.10) 生成，则

证明要点：由 \(A_kp_k=-c_k\)，\(\alpha\le1\) 时 \(\|c_k+\alpha A_kp_k\|_1=(1-\alpha)\|c_k\|_1\)，所以方向导数为 \(\nabla f_k^Tp_k-\mu^{-1}\|c_k\|_1\)。再由 (18.10) 第一行 \(\nabla f_k^Tp_k=-p_k^TW_kp_k+p_k^TA_k^T\lambda_{k+1}=-p_k^TW_kp_k-c_k^T\lambda_{k+1}\)，用 Hölder 不等式 \(|c_k^T\lambda_{k+1}|\le\|c_k\|_1\|\lambda_{k+1}\|_\infty\)。∎

所以 \(W_k\) 正定、罚权重大于最大乘子绝对值时，SQP 方向是 \(\phi_1\) 的下降方向。常用取法

与第 15 章"罚权重要超过最大乘子"的精确性阈值一致。

白话解释：\(D(\phi;p)\) 是方向导数：沿方向 \(p\) 走一小步，\(\phi\) 的变化率，即 \(\lim_{\alpha\downarrow0}[\phi(x+\alpha p)-\phi(x)]/\alpha\)。\(\phi_1\) 在 \(c_i=0\) 处有尖角、没有梯度，但方向导数总存在，所以引理用它来判断"是否下降"。(18.34) 右端两项：第一项 \(-p^TWp\) 在 \(W\) 正定时为负；第二项在 \(1/\mu>\|\lambda\|_\infty\) 时为负。

金融直觉：\(1/\mu\) 是违反约束的"罚款单价"，\(\lambda\) 是遵守约束的"影子成本"。若罚款低于遵守成本（\(1/\mu<|\lambda_i|\)），优化器会算账后宁愿交罚款、违反约束——就像监管罚金低于违规收益时，理性机构会选择违规。所以罚权重必须高于最大的影子价格，\(\ell_1\) 罚函数才"精确"，SQP 步才保证是下降方向。

18.5.2 Fletcher 价值函数

它可微。可以算出方向导数

\(\lambda_k'\) 是 \(\lambda(x)\) 的 Jacobian。约化 Hessian 正定且 \(\mu^{-1}\) 大于 (18.39) 给出的界时，\(p_k\) 是下降方向。

定理 18.3：若 \(x_k\) 不是驻点、\(Z_k^TW_kZ_k\) 正定，则 SQP 方向在 \(\mu\) 按 (18.35) 取时是 \(\phi_1\) 的下降方向，在 \(\mu\) 满足 (18.39) 时是 \(\phi_F\) 的下降方向。

罚参数更新：希望 \(\mu_k\) 最终不变。取 \(\delta>0\)：若当前 \(\mu_{k-1}^{-1}\ge\gamma+\delta\) 就保持，否则令 \(\mu_k=(\gamma+2\delta)^{-1}\)（18.40），其中 \(\ell_1\) 情形 \(\gamma=\|\lambda_{k+1}\|_\infty\)。该规则使 \(\mu\) 单调递减；若早期 \(\mu\) 变得过小，后续约束会被过度惩罚，所以有的实现加入允许 \(\mu\) 增大的启发式。

18.6 线搜索 SQP 算法

算法 18.3（非线性规划的 SQP 算法）

选 η ∈ (0, 0.5)，τ ∈ (0, 1)；选初始 (x0, λ0) 和对称正定 B0；
计算 f0, ∇f0, c0, A0；
for k = 0,1,2,...
    若满足终止测试：STOP；
    解 (18.12)（以 Bk 为 Hessian）得 pk；
    选 μk 使 pk 是 φ 在 xk 处的下降方向；
    αk = 1；
    while φ(xk + αk pk; μk) > φ(xk; μk) + η αk Dφ(xk; pk)
        αk ← τα αk，τα ∈ (0, τ)；
    x_{k+1} = xk + αk pk；
    计算 f_{k+1}, ∇f_{k+1}, c_{k+1}, A_{k+1}；
    λ_{k+1} = [A_{k+1} A_{k+1}^T]^{-1} A_{k+1} ∇f_{k+1}；（等式情形用最小二乘乘子；含不等式时取 QP 乘子）
    sk = αk pk，yk = ∇x L(x_{k+1}, λ_{k+1}) − ∇x L(xk, λ_{k+1})；
    用（阻尼）BFGS 更新 Bk 得 B_{k+1}；
end

用精确 Hessian \(W_k\) 替换 \(B_k\) 即得二阶版本（此时需注意子问题的凸性）。scipy.optimize.minimize(method="SLSQP")（Kraft 的 SLSQP 代码）基本就是这个算法：阻尼 BFGS 维护 Lagrange Hessian 近似、最小二乘形式的 QP 子问题、\(\ell_1\) 型价值函数的线搜索。

18.7 约化 Hessian SQP 方法

当二阶导难算而自由度 \(n-m\) 很小时（最优控制、化工过程控制、轨迹优化），只近似约化 Hessian 非常有效：\(M_k\) 的维数小、质量高，且即使离解较远也更可能正定。

切向收敛。保留 \(p_Y\) 时，它来自对 \(c(x)=0\) 的类牛顿迭代，通常二次收敛到零；\(p_Z\) 用拟牛顿近似，只超线性收敛。所以实践中常见

称为切向收敛。这让去掉交叉项 \(Z_k^TW_kY_kp_Y\) 变得合理；代价是收敛速度从一步超线性降为两步超线性（18.10 节），实践中差别不显著。

更新准则。由 (18.30)，\(y_k^Ts_k=p_Z^TZ_k^T\bar W_kZ_kp_Z+p_Z^TZ_k^T\bar W_kY_kp_Y\)（18.44），第一项为正，第二项符号不定；切向收敛使第一项最终占优，但不能保证，需要保护。过程 18.4（更新–跳过）：给定 \(\sum\gamma_k<\infty\) 的正数列（实践中用过 \(\gamma_k=0.1k^{-1.1}\)），只有当 \(y_k^Ts_k>0\) 且 \(\|p_Y\|\le\gamma_k\|p_Z\|\) 时才更新 \(M_k\)。跳过时说明 \(p_Y\) 相对 \(p_Z\) 不小，而 \(p_Y\) 由精确一阶信息决定，本身就能推进，\(p_Z\) 的质量此时不那么重要。

基的变化。工作集变化时约化 Hessian 维数改变：加约束时可以投影到更小的 \(M_{k+1}\)；删约束时新行列如何初始化不明显，需要若干迭代积累信息。即使维数不变，零空间基 \(Z_k\) 的选择不当也会突变——必须让 \(Z_k\) 光滑变化（18.10 节）。稳健实现相当复杂（SNOPT 等）。

算法 18.5（约化 Hessian 方法，等式约束）：每步解 \((A_kY_k)p_Y=-c_k\)、\(M_kp_Z=-Z_k^T\nabla f_k\)，\(p_k=Y_kp_Y+Z_kp_Z\)；用价值函数线搜索；在新点用 \((Y_{k+1}^TA_{k+1}^T)\lambda_{k+1}=Y_{k+1}^T\nabla f_{k+1}\) 求乘子；满足更新准则时用 (18.30) 做 BFGS 更新 \(M_k\)。

18.8 信赖域 SQP 方法

优点：能处理有效约束梯度线性相关的情形，能直接使用二阶导信息（包括不定的 Hessian）。等式情形子问题为

困难：线性化约束与信赖域可能没有交集（满足线性化约束的步都在信赖域外）。简单扩大 \(\Delta_k\) 违背信赖域的初衷；正确的观点是每步只改进可行性，仅在极限中精确满足约束。三种重构：

方法 I：约束平移（法向步 + 切向步）。先解法向子问题，忽略目标，看在信赖域的一部分内最多能多接近满足线性化约束：

解 \(v_k\) 称为法向步。再要求总步至少和法向步一样改善约束：

\(p=v_k\) 满足两个约束，所以子问题一定相容。

方法 II：两个椭球约束。\(\min\) 二次模型 s.t. \(\|A_kp+c_k\|_2\le\pi_k\)、\(\|p\|\le\Delta_k\)（18.49），\(\pi_k\) 取为某个能达到的违反量（如线性化约束残差在其 Cauchy 点处的值）。全局收敛性好，但大规模高效算法尚未建立。

方法 III：S\(\ell_1\)QP。把线性化约束以 \(\ell_1\) 罚项移进目标，只保留信赖域：

用 \(\ell_\infty\) 信赖域时，引入松弛变量即可写成标准 QP（习题 18.9）。价值函数用 \(\phi_1\) (18.52)，(18.51) 正是它的局部模型。优点：直接处理一般问题，克服不相容，放松对 \(A_k\) 正则性的要求，\(W_k\) 不必正定，全局收敛性好。缺点：对罚参数敏感，并且有 Maratos 效应。

18.8.1 实用的信赖域 SQP 算法

细化方法 I：

• 法向步用 dogleg 或 CG 近似求解。dogleg 的"牛顿点"取满足 \(A_kv+c_k=0\) 的最小范数解 \(p^B=-A_k^T(A_kA_k^T)^{-1}c_k\)，Cauchy 点沿 \(-A_k^Tc_k\)，两者都在 \(A_k^T\) 的值空间中，所以法向步 \(v_k\) 与约束切空间正交。
• 切向步：令 \(p_k=v_k+Z_ku\)，代入 (18.48) 得到切向子问题

这是标准的信赖域子问题，可用 Steihaug–CG 求解（允许 \(Z_k^TW_kZ_k\) 不定）。

• 价值函数用不平方的 \(\ell_2\)：\(\phi(x;\mu)=f(x)+\frac1\mu\|c(x)\|_2\)（18.59）；比较实际下降 ared 与预测下降 pred，按比值调整 \(\Delta_k\)。罚参数要保证预测下降至少是法向步带来的约束改善的一个比例。

算法 18.6 的主循环：算最小二乘乘子并检查 KKT 残差；解法向子问题得 \(v_k\)；算零空间基 \(Z_k\) 与 \(W_k\)；解切向子问题得 \(u_k\)；\(p_k=v_k+Z_ku_k\)；按 \(\rho_k=\mathrm{ared}/\mathrm{pred}\) 接受或拒绝步、调整 \(\Delta_k\)。scipy 的 trust-constr 方法对等式约束用的就是这种 Byrd–Omojokun 型信赖域 SQP，对不等式则套在障碍法外层（第 17 章）。

18.9 收敛速度

限于等式约束的算法 18.1。

定理 18.4：在假设 18.1 和二阶导 Lipschitz 连续的条件下，用精确 Lagrange Hessian 的局部 SQP 二次收敛。

拟牛顿版本。设投影矩阵 \(P_k=I-A_k^T(A_kA_k^T)^{-1}A_k\)（投影到约束零空间）。(18.10) 第一行乘 \(P_k\)（\(P_kA_k^T=0\)）得 \(P_kW_kp_k=-P_k\nabla f_k\)：步完全由 \(P_kW_k\) 决定，Hessian 在值空间方向的误差不影响步。

定理 18.5（Boggs–Tolle–Wang）：若拟牛顿 SQP 迭代收敛到 \(x^*\)，则它超线性收敛当且仅当

这是无约束优化 Dennis–Moré 条件（本册第 08 章）的推广：只要求 \(B_k\) 在切空间投影后沿步方向准确。

白话解释：(18.62) 的分子是"Hessian 近似误差 \(B_k-W_*\) 作用在实际步上、再投影到约束切空间后的大小"，分母是步长。比值趋于零，意思是：**不需要 \(B_k\) 整体收敛到真 Hessian，只要它在实际走的方向、且沿约束面的分量上越来越准。**垂直于约束面的方向由约束本身（\(A_kp_k=-c_k\)）决定，与 Hessian 无关，所以那里的误差无所谓。类比风险模型：你只关心组合实际持有的那些暴露方向上协方差估计准不准，没持仓的方向估计误差再大也不影响组合风险。

定理 18.6：若 \(W_*\) 正定、\(B_0\) 足够接近 \(W_*\)、\(x_0\) 足够接近 \(x^*\)，则 BFGS 版本满足 (18.62)，超线性收敛。阻尼 BFGS 只能得到较弱的 R-超线性结果。

约化 Hessian 方法：\(Z_kM_kZ_k^T\) 近似的是双侧投影 \(P_kWP_k\)，所以只能得到两步超线性收敛（定理 18.7、18.8）：\(\lim\|x_{k+2}-x^*\|/\|x_k-x^*\|=0\)。若有切向收敛，实际是一步超线性。这些结果要求零空间基光滑变化 \(\|Z_k-Z_*\|=O(\|x_k-x^*\|)\)（18.65）；任何只依赖 \(A_k\) 的 \(Z_k\) 计算过程即使 \(A_k\) 满秩也常不连续（例如 QR 中列的符号），所以实现中要让符号选择与上一步保持一致。

18.10 Maratos 效应

18.10.1 现象

价值函数可能拒绝向解推进良好的步，阻碍 SQP 的快速收敛。这由 Maratos 首先观察到。

例 18.1（Powell）：

解 \(x^*=(1,0)\)，\(\lambda^*=\frac32\)，\(\nabla^2_{xx}\mathcal L(x^*,\lambda^*)=I\)。取可行点 \(x_k=(\cos\theta,\sin\theta)\)，\(B_k=I\)。SQP 子问题的解为

\(\|x_k-x^*\|=2|\sin(\theta/2)|\)，\(\|x_k+p_k-x^*\|=2\sin^2(\theta/2)\)，比值 \(\|x_k+p_k-x^*\|/\|x_k-x^*\|^2=\frac12\)——这是一个完美的二次收敛步。但

目标和约束违反都增加了。所以任何形如 \(\phi=f+\frac1\mu h(c)\)（\(h\ge0\)，\(h(0)=0\)）的价值函数——包括 \(\ell_1\)、\(\ell_2\)、平方 \(\ell_2\)——都会拒绝这个步。原因在于：步 \(p_k\) 沿约束切线方向走，线性化约束满足了，但真实约束（圆）是弯的，\(c\) 产生了 \(O(\|p\|^2)\) 的违反；而 \(f\) 的下降量也只是 \(O(\|p\|^2)\)，罚项把它抵消掉了。不处理的话，SQP 会显著变慢，甚至失去超线性收敛。

推导拆解：几个数字怎么来的。

• \(p_k\perp x_k\)：\(p_k^Tx_k=\sin^2\theta\cos\theta-\sin\theta\cos\theta\sin\theta=0\)，所以 \(p_k\) 沿圆的切线；\(\|p_k\|=|\sin\theta|\)。
• 勾股定理：\(\|x_k+p_k\|^2=\|x_k\|^2+\|p_k\|^2=1+\sin^2\theta\)，所以 \(c(x_k+p_k)=\sin^2\theta\)，\(f(x_k+p_k)=2\sin^2\theta-(\cos\theta+\sin^2\theta)=\sin^2\theta-\cos\theta\)。
• 误差：\(x_k+p_k-x^*=(\cos\theta-\cos^2\theta,\ \sin\theta(1-\cos\theta))=(1-\cos\theta)(\cos\theta,\sin\theta)\)，长度 \(1-\cos\theta=2\sin^2(\theta/2)\)。

金融直觉：这像一个 delta 中性的对冲。沿切线走一步，一阶（delta）上约束保持满足，但约束是弯的，二阶（gamma）项让你偏离了约束面 \(\sin^2\theta\)。价值函数只看"现在偏离了多少"，不知道这个偏离是二阶小量、下一步就能修掉，于是把一个其实很好的交易否决了。二阶校正 (18.67) 相当于在 \(x_k+p_k\) 处立刻"再对冲一次 gamma"：沿法向补一个小步把约束拉回去。

18.10.2 对策

1. 用不受影响的价值函数：Fletcher 函数 \(\phi_F\) 在满足二阶充分条件的解附近会接受所有 SQP 步。
2. 二阶校正（second-order correction）：在 \(x_k+p_k\) 处重新计算约束，补一个法向步把约束违反压下去：

它是 \(A_kw+c(x_k+p_k)=0\) 的最小范数解，把 \(\|c\|\) 降到 \(O(\|x_k-x^*\|^3)\)，于是 \(x_k+p_k+w_k\) 使价值函数下降。代价是多算一次约束值。 3. 非单调策略（watchdog）：允许价值函数偶尔上升（"松弛步"），若在若干步内没有获得充分下降，就回到松弛步之前的点做正常的线搜索。

**算法 18.7（带二阶校正的 SQP）**的逻辑：先试 \(\alpha=1\)；若 \(\phi_1\) 不满足充分下降，计算 \(w_k\) 并试 \(x_k+p_k+w_k\)；成功则接受，否则回溯缩短 \(\alpha\)。可以证明有限步后 \(\alpha_k=1\) 总能产生新迭代点，价值函数不再干扰，获得与局部 SQP 相同的超线性收敛。

**算法 18.8（watchdog）**的逻辑：若完整 SQP 步满足充分下降就正常接受；否则暂时接受它（松弛步），再走一步线搜索 SQP；若两步后的点相对原点有充分下降就继续，否则退回原点做线搜索。至少三分之一的迭代获得充分下降，可证局部超线性收敛。实践中常允许 5–10 步上升。相对二阶校正，它可能节省约束求值，两者优劣尚无定论。

18.11 量化实战

18.11.1 在哪里用

• 非线性约束的组合优化：波动率目标 \(\sqrt{w^T\Sigma w}\le\sigma^*\)、风险贡献上限 \(w_i(\Sigma w)_i\le c\,w^T\Sigma w\)、市场冲击成本 \(\sum\eta_i|w_i-w_i^0|^{3/2}\)、CVaR/下行风险约束——问题不再是 QP。中小规模直接用 SLSQP（本章算法 18.3）；需要乘子或大规模时用 trust-constr（18.8.1 节 + 第 17 章障碍法）、IPOPT、SNOPT、KNITRO。
• 读懂 SLSQP 的报错：status 4 "Inequality constraints incompatible" 是 18.3.3 节的线性化不相容；status 8 "Positive directional derivative for linesearch" 通常是 Hessian 近似不好或价值函数拒绝步（含 Maratos 效应）；status 9 "Iteration limit reached" 在不可行问题上很常见。对策：检查可行性、换起点、做变量缩放、软化可能冲突的约束。
• 模型校准：Heston、SABR 校准带非线性约束（如 Feller 条件 \(2\kappa\theta>\sigma^2\)）时可用 SQP；精确 Hessian 往往拿不到，阻尼 BFGS 是标准选择。
• 最优执行：Almgren–Chriss 型执行轨迹加入非线性冲击或风险约束后是结构化的 NLP，自由度少，适合约化 Hessian SQP 或利用带状结构的内点法。
• 乘子 = 影子价格：SQP 返回的 QP 乘子告诉你放宽某个风险约束能换来多少收益。
• 热启动：日频再平衡时以前一日的解作初值，SQP 通常几步收敛。

18.11.2 代码一：局部 SQP 的二次收敛与 Maratos 效应

用原书习题 18.2 的问题实现算法 18.1（精确 Hessian）：

import numpy as np

# ---------- (a) 算法 18.1 局部 SQP（精确 Hessian），原书习题 18.2 ----------
def prod_except(x, idx):
    return np.prod(np.delete(x, idx))
def f_all(x):
    P = np.prod(x); e = np.exp(P); g = x[0]**3 + x[1]**3 + 1
    gP = np.array([prod_except(x, [i]) for i in range(5)])
    HP = np.array([[0 if i == j else prod_except(x, [i, j]) for j in range(5)] for i in range(5)])
    dg = np.array([3*x[0]**2, 3*x[1]**2, 0, 0, 0]); Hg = np.diag([6*x[0], 6*x[1], 0, 0, 0])
    f = e - 0.5 * g**2
    grad = e * gP - g * dg
    hess = e * (np.outer(gP, gP) + HP) - (np.outer(dg, dg) + g * Hg)
    return f, grad, hess
def c_all(x):
    c = np.array([x @ x - 10, x[1]*x[2] - 5*x[3]*x[4], x[0]**3 + x[1]**3 + 1])
    A = np.array([2 * x,
                  [0, x[2], x[1], -5*x[4], -5*x[3]],
                  [3*x[0]**2, 3*x[1]**2, 0, 0, 0]])
    H1 = 2 * np.eye(5)
    H2 = np.zeros((5, 5)); H2[1, 2] = H2[2, 1] = 1; H2[3, 4] = H2[4, 3] = -5
    H3 = np.diag([6*x[0], 6*x[1], 0, 0, 0])
    return c, A, [H1, H2, H3]

x = np.array([-1.71, 1.59, 1.82, -0.763, -0.763]); lam = np.zeros(3)
print(" k   ||∇L||          ||c||           ||Δx||")
for k in range(10):
    f, g, Hf = f_all(x); c, A, Hc = c_all(x)
    W = Hf - sum(l * H for l, H in zip(lam, Hc))          # 拉格朗日 Hessian
    KKT = np.block([[W, -A.T], [A, np.zeros((3, 3))]])
    sol = np.linalg.solve(KKT, np.r_[-g, -c])             # (18.10)
    p, lam_new = sol[:5], sol[5:]
    print(f"{k:2d}  {np.linalg.norm(g - A.T @ lam):.3e}   {np.linalg.norm(c):.3e}   {np.linalg.norm(p):.3e}")
    if np.linalg.norm(p) < 1e-14: break
    x, lam = x + p, lam_new
print("x* =", x.round(6), " λ* =", lam.round(6))

# ---------- (b) Maratos 效应：例 18.1 ----------
print("\n== Maratos 效应（例 18.1）==")
fM = lambda x: 2 * (x @ x - 1) - x[0]
cM = lambda x: x @ x - 1
for theta in [0.5, 0.1, 0.01]:
    xk = np.array([np.cos(theta), np.sin(theta)])
    p = np.array([np.sin(theta)**2, -np.sin(theta) * np.cos(theta)])   # (18.66)
    xs = np.array([1.0, 0.0]); x1 = xk + p
    A = 2 * xk[None, :]
    w = -A.T @ np.linalg.solve(A @ A.T, [cM(x1)])                     # 二阶校正 (18.67)
    x2 = x1 + w
    phi = lambda z, inv_mu=2.0: fM(z) + inv_mu * abs(cM(z))           # 1/μ = 2 > |λ*| = 1.5
    print(f"θ={theta:4.2f}: 误差 {np.linalg.norm(xk-xs):.2e} → SQP步 {np.linalg.norm(x1-xs):.2e} "
          f"→ 校正后 {np.linalg.norm(x2-xs):.2e} | φ1: {phi(xk):.6f} → {phi(x1):.6f}"
          f"（{'上升' if phi(x1) > phi(xk) else '下降'}） → {phi(x2):.6f}"
          f"（{'上升' if phi(x2) > phi(xk) else '下降'}）")

输出：

 k   ||∇L||          ||c||           ||Δx||
 0  4.069e-01   7.563e-02   1.184e-02
 1  2.871e-03   1.763e-04   1.721e-03
 2  1.104e-05   3.249e-06   2.729e-06
 3  1.128e-11   8.046e-12   1.331e-11
 4  8.327e-17   4.441e-16   5.950e-16
x* = [-1.717144  1.59571   1.827246 -0.763643 -0.763643]  λ* = [-0.040163  0.037958 -0.005223]

== Maratos 效应（例 18.1）==
θ=0.50: 误差 4.95e-01 → SQP步 1.22e-01 → 校正后 7.49e-03 | φ1: -0.877583 → -0.188036（上升） → -0.953745（下降）
θ=0.10: 误差 1.00e-01 → SQP步 5.00e-03 → 校正后 1.25e-05 | φ1: -0.995004 → -0.965104（上升） → -0.999913（下降）
θ=0.01: 误差 1.00e-02 → SQP步 5.00e-05 → 校正后 1.25e-09 | φ1: -0.999950 → -0.999650（上升） → -1.000000（下降）

解读：

1. 二次收敛：KKT 残差 \(4\times10^{-1}\to3\times10^{-3}\to10^{-5}\to10^{-11}\to10^{-16}\)，有效数字每步翻倍；即使初始乘子取 \(\lambda_0=0\)，4 步就达到机器精度。解 \(x^*\approx(-1.717,1.596,1.827,-0.764,-0.764)\)，与原书给出的近似解一致。
2. Maratos 效应：\(\theta=0.1\) 时 SQP 步把误差从 \(10^{-1}\) 降到 \(5\times10^{-3}\)（\(\approx\frac12\|x_k-x^*\|^2\)，二次收敛），但 \(\ell_1\) 价值函数从 \(-0.995\) 升到 \(-0.965\)，线搜索会拒绝这个好步。罚参数 \(1/\mu=2\) 已经超过 \(|\lambda^*|=1.5\)，所以这不是罚参数太小的问题。
3. 二阶校正：加上 \(w_k\) 后误差进一步降到 \(1.25\times10^{-5}\)，价值函数下降，步被接受。代价只是在 \(x_k+p_k\) 处多算一次约束。

18.11.3 代码二：波动率目标 + 风险贡献上限的组合（SLSQP 与 trust-constr）

问题：25 只资产，最大化预期收益 \(\alpha^Tw\)，约束为预算 \(\mathbf 1^Tw=1\)、长仓、个股上限 15%、组合波动 ≤ 12%、单资产风险贡献占比 ≤ 10%：\(w_i(\Sigma w)_i\le0.1\,w^T\Sigma w\)。后两条是非线性约束（二次，且风险贡献约束非凸）。

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

rng = np.random.default_rng(8)
n = 25
B = rng.standard_normal((n, 3)) * np.array([0.15, 0.08, 0.06]) + np.array([0.18, 0, 0])
Sigma = B @ B.T + np.diag(rng.uniform(0.01, 0.06, n))
alpha = 0.03 + 0.04 * rng.standard_normal(n)
sig_tgt, rc_cap, ub = 0.12, 0.10, 0.15      # 波动目标 12%，单资产风险贡献 ≤10%，权重 ≤15%

f    = lambda w: -alpha @ w
gf   = lambda w: -alpha
vol2 = lambda w: w @ Sigma @ w
def rc_cons(w):                       # rc_cap * w'Σw - w_i (Σw)_i >= 0
    s = Sigma @ w
    return rc_cap * (w @ s) - w * s
def rc_jac(w):
    s = Sigma @ w
    return rc_cap * 2 * s[None, :] - (np.diag(s) + w[:, None] * Sigma)
cons = [{"type": "eq",   "fun": lambda w: w.sum() - 1, "jac": lambda w: np.ones((1, n))},
        {"type": "ineq", "fun": lambda w: sig_tgt**2 - vol2(w), "jac": lambda w: -2 * Sigma @ w},
        {"type": "ineq", "fun": rc_cons, "jac": rc_jac}]
w0 = np.full(n, 1 / n)
r = minimize(f, w0, jac=gf, method="SLSQP", bounds=[(0, ub)] * n, constraints=cons,
             options={"ftol": 1e-12, "maxiter": 500})
w = r.x
print(f"SLSQP: {r.message} 迭代 {r.nit}，预期收益 {alpha @ w:.4%}，波动 {np.sqrt(vol2(w)):.4%}")
print(f"  最大风险贡献占比 {np.max(w * (Sigma @ w)) / vol2(w):.4f}，持仓数 {(w > 1e-6).sum()}，"
      f"顶格 {(w > ub - 1e-6).sum()}")

# trust-constr：同一问题，返回拉格朗日乘子（影子价格）
r2 = minimize(f, w0, jac=gf, hess=lambda w: np.zeros((n, n)), method="trust-constr", bounds=Bounds(0, ub),
              constraints=[LinearConstraint(np.ones((1, n)), 1, 1),
                           NonlinearConstraint(vol2, -np.inf, sig_tgt**2, jac=lambda w: 2 * Sigma @ w,
                                               hess=lambda w, v: v[0] * 2 * Sigma),
                           NonlinearConstraint(rc_cons, 0, np.inf, jac=rc_jac)],
              options={"gtol": 1e-10, "xtol": 1e-12, "maxiter": 5000})
print(f"trust-constr: 预期收益 {alpha @ r2.x:.4%}，与 SLSQP 权重最大差 {np.abs(r2.x - w).max():.1e}，迭代 {r2.nit}，{r2.message}")
lam_vol = abs(r2.v[1][0])
print(f"  方差约束乘子 |λ| = {lam_vol:.4f}  ⇒ dR/dσ = 2σλ = {2*sig_tgt*lam_vol:.4f}："
      f"波动上限每放宽 1 个百分点，预期收益约增 {2*sig_tgt*lam_vol*0.01:.3%}")

# 有限差分验证影子价格：把方差上限放宽 1e-4
def solve(s2):
    c2 = [cons[0], {"type": "ineq", "fun": lambda w: s2 - vol2(w), "jac": lambda w: -2 * Sigma @ w}, cons[2]]
    return minimize(f, w, jac=gf, method="SLSQP", bounds=[(0, ub)] * n, constraints=c2,
                    options={"ftol": 1e-14, "maxiter": 500})
d = 1e-4
print(f"  有限差分 d(收益)/d(方差上限) = {(-solve(sig_tgt**2 + d).fun + solve(sig_tgt**2).fun) / d:.4f}")

# 不可行问题：要求波动 ≤ 2%（低于最小方差组合的波动）
bad = [cons[0], {"type": "ineq", "fun": lambda w: 0.02**2 - vol2(w), "jac": lambda w: -2 * Sigma @ w}]
mv = minimize(vol2, w0, jac=lambda w: 2 * Sigma @ w, method="SLSQP", bounds=[(0, ub)] * n,
        constraints=[cons[0]], options={"ftol": 1e-14})
print(f"\n最小方差组合的波动 = {np.sqrt(mv.fun):.2%}")
r3 = minimize(f, w0, jac=gf, method="SLSQP", bounds=[(0, ub)] * n, constraints=bad)
print("波动上限 2% 时 SLSQP 返回:", r3.status, r3.message)

输出：

SLSQP: Optimization terminated successfully 迭代 15，预期收益 5.3006%，波动 12.0000%
  最大风险贡献占比 0.1000，持仓数 13，顶格 3
trust-constr: 预期收益 5.2999%，与 SLSQP 权重最大差 3.1e-04，迭代 56，`gtol` termination condition is satisfied.
  方差约束乘子 |λ| = 1.0388  ⇒ dR/dσ = 2σλ = 0.2493：波动上限每放宽 1 个百分点，预期收益约增 0.249%
  有限差分 d(收益)/d(方差上限) = 1.0312

最小方差组合的波动 = 5.83%
波动上限 2% 时 SLSQP 返回: 9 Iteration limit reached

解读：

1. SLSQP 15 次迭代求得最优组合：波动恰好 12%，最大风险贡献恰好 10%——两条非线性约束都是有效约束；13 只持仓中 3 只顶格。
2. trust-constr 给出几乎相同的解（权重差 \(3\times10^{-4}\)，来自其障碍参数的终止精度），并返回乘子。方差约束的乘子 1.039 与 SLSQP 有限差分 1.031 基本一致（差异来自 trust-constr 的终止精度和有限差分步长）。换算成波动率：\(dR/d\sigma=2\sigma\lambda\approx0.25\)，即波动上限每放宽 1 个百分点，预期收益约增 25bp——这是"风险预算的边际收益"，可以直接用于风险预算在各策略间的分配。

推导拆解：\(dR/d\sigma=2\sigma\lambda\) 是链式法则。约束写的是方差上限 \(v=\sigma^2\)，乘子 \(\lambda\) 是最优收益对 \(v\) 的边际：\(dR/dv=\lambda\)（包络定理，乘子 = 影子价格）。又 \(dv/d\sigma=2\sigma\)，所以 \(dR/d\sigma=\frac{dR}{dv}\cdot\frac{dv}{d\sigma}=2\sigma\lambda=2\times0.12\times1.039\approx0.249\)。代码取 abs(r2.v[1][0]) 是因为 trust-constr 的乘子符号约定与本书 \(\mathcal L=f-\lambda^Tc\) 不同，这里只关心大小。同理也可以换算成"收益/风险"的边际斜率，用于在多个策略间分配风险预算：边际斜率高的策略应多分风险。

3. 不可行问题的表现：波动上限 2% 低于最小方差组合的 5.83%，问题不可行。SLSQP 没有弹性模式，只是跑满迭代次数（status 9），并不明确告诉你"不可行"。实务中应先用最小方差等辅助问题检查可行性，再做正式优化。

本章小结

SQP 用二次子问题逼近约束问题：等式情形下它等价于对 KKT 方程组的牛顿法，新 \(x\) 是 QP 的解、新 \(\lambda\) 是 QP 的乘子，局部二次收敛；不等式情形线性化全部约束（IQP）或先选工作集（EQP），接近解时有效集稳定。步计算可用增广系统、值空间或零空间法；最小二乘乘子把迭代变为纯原始；约化 Hessian 方法去掉交叉项，依靠切向收敛只损失到两步超线性。Lagrange Hessian 常不正定，阻尼 BFGS 用 \(r_k=\theta y_k+(1-\theta)B_ks_k\) 保证正定（SLSQP 即如此），信赖域法可用 SR1 或精确 Hessian，增广拉格朗日 Hessian 和约化 Hessian 是另两种凸化途径。线搜索 SQP 用 \(\ell_1\) 价值函数时罚权重须超过 \(\|\lambda\|_\infty\)；信赖域 SQP 用法向步 + 切向步、双椭球或 S\(\ell_1\)QP 解决线性化约束与信赖域不相容的问题。拟牛顿 SQP 超线性收敛的充要条件是投影的 Dennis–Moré 条件。Maratos 效应使 \(f+h(c)/\mu\) 型价值函数拒绝二次收敛的好步，对策是 Fletcher 函数、二阶校正 \(w_k=-A^T(AA^T)^{-1}c(x_k+p_k)\) 或 watchdog。

练习

基础

1. 证明定理 18.4（用牛顿法的二次收敛定理）。
2. 证明阻尼 BFGS 在 \(\theta_k\ne1\) 时满足 \(s_k^Tr_k=0.2\,s_k^TB_ks_k\)，从而 \(B_{k+1}\) 正定。
3. 写出约束 \(x_1^2+x_2^2=1\) 在 \((0,0)\)、\((0,1)\)、\((0.1,0.02)\)、\(-(0.1,0.02)\) 处的线性化，指出哪些点线性化退化或与信赖域 \(\|p\|\le0.5\) 不相容。
4. 对例 18.1 验证 (18.66)：写出 SQP 子问题并求解。
5. 证明 (18.51) 在 \(\ell_\infty\) 信赖域下可以写成一个 QP。 提示：对每个 \(|\cdot|\) 引入 \(u_i-v_i\) 拆分，对 \([\cdot]^-\) 引入非负变量 \(t_i\ge-(\cdot)\)。

进阶

6. 证明 (18.38)。
7. 实现约化 Hessian SQP（算法 18.5）求解习题 18.2 的问题，比较它与 18.11.2 中精确 Hessian SQP 的迭代次数，观察 \(\|p_Y\|/\|p_Z\|\) 是否趋于零（切向收敛）。
8. 在 18.11.2 的 Maratos 例子上实现完整的线搜索 SQP（\(\ell_1\) 价值函数），分别用 (a) 纯回溯、(b) 二阶校正（算法 18.7）、(c) watchdog（\(\hat t=1\)）从 \(\theta=1\) 出发，比较迭代次数和收敛速度。
9. 用 CG 解切向子问题时，证明 \(Z_k\) 只以 \(Z_k(Z_k^TZ_k)^{-1}Z_k^Tw\) 的形式出现，它等于 \(w-A_k^T(A_kA_k^T)^{-1}A_kw\)，因此无需显式构造零空间基。
10. 在 18.11.3 的问题中把风险贡献上限从 10% 逐步收紧到 5%，记录 SLSQP 的迭代次数、返回状态和解；当 \(1/n=4\%\) 附近时约束接近要求"等风险贡献"，观察问题难度如何变化，并与第 17 章的障碍型风险平价凸化比较。

原书推荐习题：18.2（实现局部 SQP，体会二次收敛）、18.3（阻尼 BFGS 正定性）、18.6（线性化约束不相容的直观认识）、18.9（S\(\ell_1\)QP 子问题的 QP 化）、18.10（不显式构造零空间基的投影 CG）；并建议手算例 18.1，亲自验证 Maratos 效应与二阶校正。原书第 18 章共 18.1–18.10 题。

原书对照

更正说明：算法 18.3 与 18.5 中乘子公式原文多一个负号，按 (18.15)/(18.16) 应无负号；(18.61) 分子按推导应为切向模型的增加量；图 18.4 取 \(\theta=\pi/2\) 时起点应为 \((0,1)\)；习题 18.10 中零空间投影应为 \(w-A_k^T(A_kA_k^T)^{-1}A_kw\)。原书正文页码 = PDF 页码减 18（已用 PDF 页眉核对，第 17 章至附录均如此）。