第 14 章 线性规划：内点法

学习目标

读完本章，你应当能够：

1. 把 LP 的 KKT 条件写成非线性方程组 \(F(x,\lambda,s)=0\) 加非负约束，并说明为什么必须让迭代点严格保持 \((x,s)>0\)。
2. 理解中心路径、对偶度量 \(\mu=x^Ts/n\) 与中心化参数 \(\sigma\) 的含义，说清仿射尺度方向（\(\sigma=0\)）与中心化方向（\(\sigma=1\)）的取舍。
3. 写出路径跟踪方法的两种邻域 \(\mathcal N_2(\theta)\)、\(\mathcal N_{-\infty}(\gamma)\)，读懂长步路径跟踪算法 \(O(n\log1/\epsilon)\) 复杂度证明的主线。
4. 按步骤实现 Mehrotra 预测–校正算法，并知道它是绝大多数商业内点求解器的基础。
5. 理解每步线性代数的两种形式（增广系统、正规方程 \(AD^2A^T\)），知道稠密列和后期病态带来的问题。
6. 能读懂内点求解器日志（\(\mu\)、原始/对偶残差、对偶间隙），并在实际 LP 中比较单纯形法与内点法。

读前导读

这一章在解决什么问题

上一章的单纯形法沿着可行区域的"棱角"一个顶点一个顶点地走，问题一大，要走的顶点就成千上万。本章换一个思路：不走边界，从可行区域内部直接穿过去，逼近最优点。做法是把 LP 的 KKT 条件（原始可行、对偶可行、互补松弛）看成一组方程，用牛顿法求解，同时始终保持变量严格为正。

牛顿法你其实很熟。金融计算器求 IRR 或到期收益率时，就是在反复"用切线近似、解出下一个猜测值"，几步就收敛，这就是一元牛顿法。本章把它推广到几百上千个方程同时求解。难点不在牛顿法本身，而在"互补松弛 \(x_is_i=0\)"这类条件：直接解会掉进错误的解（变量变负）。内点法的对策是先把目标放宽成"\(x_is_i\) 都等于一个小正数 \(\tau\)"，再让 \(\tau\) 逐步降到零。这条随 \(\tau\) 移动的轨迹叫中心路径。

可以把它想成风控里的"软限额"：不是到了限额才硬性禁止，而是越接近限额惩罚越重，于是组合自然停在限额内部；再逐步减轻惩罚，组合就慢慢贴近真正的最优边界。读完本章你还能看懂 MOSEK、Gurobi 等求解器的 barrier 日志：每行打印的 \(\mu\) 就是对偶间隙（离最优还有多远）的平均值。

需要先想起来的数学

1. 牛顿法解方程组。 一元情形：解 \(g(x)=0\)，在当前点做切线 \(g(x)+g'(x)\Delta x=0\)，得 \(\Delta x=-g(x)/g'(x)\)。例：解 \(x^2-2=0\)，从 \(x=1\) 出发，\(\Delta x=-(1-2)/2=0.5\)，得 1.5，再一步得 1.4167，已很接近 \(\sqrt2\)。多元情形把 \(g'\) 换成 Jacobian 矩阵 \(J\)（第 \(i\) 个方程对第 \(j\) 个变量的偏导排成的矩阵），解线性方程组 \(J\Delta=-F\)。见 第 00 册第 05 章 多元微积分与优化。

2. 对角矩阵记号。 \(X=\mathrm{diag}(x)\) 是把向量 \(x\) 放在对角线上、其余为 0 的方阵；\(e\) 是全 1 向量。于是 \(XSe\) 就是逐分量乘积组成的向量 \((x_1s_1,\dots,x_ns_n)^T\)，\(X^{-1}\) 就是对角线取倒数。这只是一种紧凑写法，没有新内容。见 第 00 册第 06 章 线性代数速成。

3. 范数。 \(\|v\|_2=\sqrt{\sum v_i^2}\) 是通常的长度；\(\|v\|_1=\sum|v_i|\)；\(\|v\|_\infty=\max|v_i|\)。例：\(v=(3,-4)\)，三者分别是 5、7、4。本章用范数衡量"离中心路径有多远"。

4. 对数与几何衰减。 若每步把误差乘以 \((1-\delta/n)\)，\(k\) 步后误差约为原来的 \(e^{-k\delta/n}\)。要从 \(\mu_0\) 降到 \(\epsilon\)，需要 \(k\approx\frac n\delta\ln\frac{\mu_0}\epsilon\) 步。这和"按固定比例摊销的贷款余额多久降到某个水平"是同一种计算。记号 \(O(\cdot)\) 表示"同阶增长，忽略常数"。见 第 00 册第 04 章 级数与收敛 和 第 07 章 概率中的分析工具。

怎么读这一章

核心必读是 14.2（原始–对偶框架：KKT 方程组、牛顿步、中心路径、\(\mu\) 和 \(\sigma\)）、14.4（Mehrotra 预测–校正，实际软件用的就是它）和 14.7（代码与求解器日志）。建议先读 14.2，再直接跳到 14.4 和 14.7.2，对照代码输出理解每一步。

第一次可以跳过的是 14.3.1 的复杂度证明（只需记住结论：迭代次数随问题规模增长很慢）、14.3.2 的其他方法简介和 14.6 的推广。14.5 的线性代数如果你不打算自己写求解器，看懂"每步的主要成本是解一个 \(AD^2A^T\) 方程组，稠密的约束行会让它变慢"这一点就够了。

14.1 动机：从内部逼近解

先说结论：原始–对偶内点法（primal–dual interior-point methods）对 LP 的 KKT 条件施用牛顿法，同时让所有迭代点严格位于非负象限内部。它每步代价比单纯形法高，但迭代次数很少（通常 20–60 次），且几乎不随问题规模增长，在大规模问题上与单纯形法强有力地竞争。

历史线索。单纯形法最坏情况下是指数时间（第 13 章 Klee–Minty 例子）。Khachiyan 的椭球法在最坏情况下是多项式的，但它在几乎所有问题上都接近最坏界，毫无竞争力。1984 年 Karmarkar 提出的投影算法既是多项式的，又有不错的实际表现，引发了大量研究：仿射尺度、对数障碍、势函数下降、路径跟踪、原始–对偶、不可行内点等方法相继出现。到 90 年代初，原始–对偶方法脱颖而出，成为最有效的实用方法。

两类方法的几何区别：单纯形法沿可行多面体的边界逐个检查顶点；内点法只在极限时才接近边界，可以从内部甚至外部（不可行内点法）逼近解，但从不落在边界上。

14.2 原始–对偶框架

14.2.1 KKT 条件写成方程组

原始问题与对偶问题：

本章对偶变量记为 \(\lambda\)（第 13 章记为 \(\pi\)）。KKT 条件：

写成映射 \(F:\mathbb R^{2n+m}\to\mathbb R^{2n+m}\)：

其中 \(X=\mathrm{diag}(x)\)，\(S=\mathrm{diag}(s)\)，\(e=(1,\dots,1)^T\)。

三个等式中前两个线性、第三个只是"轻度非线性"（双线性），单独求解并不难。非负性要求才是内点法一切复杂性的来源。如果只解 \(F=0\) 而不管符号，会碰到大量伪解（spurious solutions）：满足 \(F=0\) 但 \((x,s)\) 有负分量的点，它们与 LP 的解毫无关系。

例（习题 14.1）：\(\min x_1\) s.t. \(x_1+x_2=1\)，\(x\ge0\)。真正的原始–对偶解是 \(x^*=(0,1)\)，\(\lambda^*=0\)，\(s^*=(1,0)\)。但 \(x=(1,0)\)，\(\lambda=1\)，\(s=(0,-1)\) 也满足 \(F=0\)：\(A^T\lambda+s=(1,1)+(0,-1)=(1,0)=c\)，\(x_1+x_2=1\)，\(x_1s_1=x_2s_2=0\)。它只是违反了 \(s\ge0\)。

所以内点法要求迭代点严格满足 \(x>0,\ s>0\)，"内点"之名即由此而来。定义可行集与严格可行集：

14.2.2 牛顿步与仿射尺度方向

对 \(F=0\) 用牛顿法：\(J(x,\lambda,s)\,(\Delta x,\Delta\lambda,\Delta s)=-F\)。若当前点严格可行（前两块残差为零），牛顿方程为

推导拆解：(14.7) 的系数矩阵就是 \(F\) 的 Jacobian，逐块求导即可。第一块 \(A^T\lambda+s-c\)：对 \(x\) 求导为 0，对 \(\lambda\) 为 \(A^T\)，对 \(s\) 为单位阵 \(I\)。第二块 \(Ax-b\)：只对 \(x\) 有导数 \(A\)。第三块的第 \(i\) 个分量 \(x_is_i\)：对 \(x_i\) 求导得 \(s_i\)，对 \(s_i\) 求导得 \(x_i\)，排成矩阵就是 \(S\) 和 \(X\)。 第三行的含义最直观：\((x_i+\Delta x_i)(s_i+\Delta s_i)=x_is_i+s_i\Delta x_i+x_i\Delta s_i+\Delta x_i\Delta s_i\)。牛顿法丢掉最后的二阶项，令前三项等于目标值 0，就得到 \(s_i\Delta x_i+x_i\Delta s_i=-x_is_i\)。被丢掉的 \(\Delta x_i\Delta s_i\) 后面会在 Mehrotra 校正步里被补回来。

这个纯牛顿方向称为仿射尺度方向（affine scaling direction）。问题在于：沿它走全步通常会违反 \((x,s)\ge0\)，只能走很小的步 \(\alpha\ll1\)，进展有限。原始–对偶方法做两项修改：

1. 让方向偏向非负象限内部，以便走得更远；
2. 防止 \((x,s)\) 的分量过早靠近边界。

14.2.3 中心路径

中心路径（central path） \(\mathcal C\) 是由参数 \(\tau>0\) 刻画的一条严格可行点曲线，每点 \((x_\tau,\lambda_\tau,s_\tau)\) 满足

它与 KKT 条件唯一的差别是互补条件的右端从 \(0\) 变成 \(\tau\)：要求所有乘积 \(x_is_i\) 相等。可以证明，\(\mathcal F^o\) 非空时每个 \(\tau>0\) 唯一确定中心路径上的一点；\(\tau\to0\) 时 (14.8) 越来越接近 (14.3)，若 \(\mathcal C\) 收敛，必收敛到 LP 的原始–对偶解。

直观上，中心路径是一条"避开伪解的安全通道"：沿着它走，\(x,s\) 始终严格为正，所有 \(x_is_i\) 以大致相同的速率降到零，不会有某个分量过早撞到边界。

白话解释：回忆第 13 章 (13.6)：任一可行方案的成本 = 对偶给出的资源总价值 + "浪费项" \(x^Ts=\sum_ix_is_i\)。最优就是浪费为零。KKT 的互补条件要求每一项 \(x_is_i\) 都为零，但"为零"有两种方式（\(x_i=0\) 或 \(s_i=0\)），直接解方程分不清哪种是对的，还可能用负数凑出零，这就是上面的伪解。 中心路径的做法是暂时不要求为零，而是要求每一项都等于同一个正数 \(\tau\)。这时 \(x_i\) 和 \(s_i\) 都必须严格为正，伪解自然被排除。然后让 \(\tau\) 从大到小逐步降到 0：\(x_i\) 和 \(s_i\) 中"该为零的那个"会慢慢趋于零，"不该为零的那个"留在正值。哪个该为零，交给算法在过程中自然决定，不需要像单纯形法那样一个个去猜。 金融上的类比是软限额：与其在限额处设一堵硬墙，不如设一个越靠近越贵的惩罚（对数障碍 \(-\tau\log x_i\)，\(x_i\to0\) 时趋于无穷），组合就会停在墙内；惩罚系数 \(\tau\) 逐步减小，组合就逐步贴近真正的最优边界。

第 17 章将看到，中心路径正是对数障碍函数 \(c^Tx-\tau\sum\log x_i\) 的极小点轨迹；\(x_is_i=\tau\) 就是障碍问题的一阶条件。

14.2.4 对偶度量与中心化参数

定义对偶度量（duality measure）

对可行点，\(x^Ts=c^Tx-b^T\lambda\) 就是对偶间隙，所以 \(\mu\) 直接衡量离最优还有多远。原始–对偶方法不朝 \(F=0\) 走纯牛顿步，而是朝中心路径上 \(\tau=\sigma\mu\) 的点走牛顿步，\(\sigma\in[0,1]\) 称为中心化参数（centering parameter）：

• \(\sigma=1\)：中心化方向，朝 \(x_is_i=\mu\) 的中心点走，强烈偏向内部，几乎不降低 \(\mu\)，但为下一步走大步打下基础；
• \(\sigma=0\)：标准牛顿步，即仿射尺度方向；
• 多数算法取 \(\sigma\in(0,1)\)，在降低 \(\mu\) 与改善中心性之间权衡。

白话解释：\(\mu\) 是平均每个变量的"浪费"，\(n\mu\) 就是对偶间隙：当前方案的成本与对偶下界之差，相当于"最多还能再省多少"。\(\sigma\) 决定这一步的目标是把 \(\mu\) 降到多少：\(\sigma=0\) 是一步到位地追求零间隙，\(\sigma=1\) 是间隙不变、只把各个 \(x_is_i\) 拉平。 为什么不总是取 \(\sigma=0\)？因为如果某些 \(x_is_i\) 远小于平均值，它们离边界太近，沿牛顿方向稍走一步就会变负，只能走极小的步，反而进展很慢。适度中心化（\(\sigma>0\)）先把落后的分量拉回来，下一步才能放心走大步。这有点像再平衡：先把偏离太大的头寸拉回目标附近，后续调整才有空间。

Framework 14.1（原始–对偶框架）：给定 \((x^0,\lambda^0,s^0)\in\mathcal F^o\)。对 \(k=0,1,\dots\)：选 \(\sigma_k\in[0,1]\)，以 \(\mu_k=(x^k)^Ts^k/n\) 解 (14.11) 得方向；取步长 \(\alpha_k\) 使 \((x^{k+1},s^{k+1})>0\)，更新。\(\sigma_k\) 与 \(\alpha_k\) 的选法决定了具体算法及其理论性质。

14.2.5 不可行内点法

严格可行的初始点通常很难找。**不可行内点法（infeasible-interior-point methods）**只要求 \(x^0,s^0>0\)。定义残差

步方程改为

由于前两块方程是线性的，若某步取全步 \(\alpha=1\)，残差立刻归零，此后迭代保持可行；取 \(\alpha<1\) 时残差按 \((1-\alpha)\) 的比例缩小。实用软件都用这种形式。

14.3 路径跟踪方法

**路径跟踪方法（path-following methods）**显式地把迭代点限制在中心路径的某个邻域内，沿 \(\mathcal C\) 走向解，以 \(\mu\) 作为进展的度量，迫使 \(\mu_k\to0\)。两种常用邻域：

典型取值 \(\theta=0.5\)、\(\gamma=10^{-3}\)。\(\mathcal N_{-\infty}(\gamma)\) 只要求每个乘积不低于平均值的 \(\gamma\) 倍，非常宽松，\(\gamma\) 很小时几乎覆盖整个可行域；\(\mathcal N_2(\theta)\) 严格得多。与一般非线性方程的同伦法相比，同伦法停留在路径的一个"细管"里，而原始–对偶方法的邻域是锥形的：\(\mu\) 大时宽松，\(\mu\to0\) 时收窄。

白话解释：两个邻域都在问"各个 \(x_is_i\) 是否足够均匀"。\(\mathcal N_2(\theta)\) 要求它们与平均值 \(\mu\) 的整体偏差（欧氏距离）不超过 \(\theta\mu\)；\(\mathcal N_{-\infty}(\gamma)\) 只要求没有哪一项低于平均值的 \(\gamma\) 倍，相当于"只设下限、不设上限"的集中度规则。例：\(n=2\)、\(x_1s_1=0.1\)、\(x_2s_2=1.9\)，则 \(\mu=1\)，\(\|XSe-\mu e\|_2=\sqrt{0.81+0.81}\approx1.27\)，不在 \(\mathcal N_2(0.5)\) 里；但 \(0.1\ge10^{-3}\times1\)，在 \(\mathcal N_{-\infty}(10^{-3})\) 里。同伦法是解非线性方程的一类方法：从一个容易解的方程出发，逐步变形到目标方程，沿途跟踪解的轨迹。

Algorithm 14.2（长步路径跟踪，long-step path-following）：给定 \(\gamma\in(0,1)\)，\(0<\sigma_{\min}<\sigma_{\max}<1\)，初始点 \(\in\mathcal N_{-\infty}(\gamma)\)。每步选 \(\sigma_k\in[\sigma_{\min},\sigma_{\max}]\)，解 (14.11)，取 \(\alpha_k\) 为 \([0,1]\) 中使新点仍在 \(\mathcal N_{-\infty}(\gamma)\) 内的最大值，更新。

\(\sigma_{\min}>0\) 保证每个方向一开始先改善中心性（离开邻域边界），所以总能走一个不太小的步。下面证明它是多项式算法。

14.3.1 复杂度分析

整条证明链只用初等不等式：引理 14.1 → 引理 14.2 → 定理 14.3 → 定理 14.4。

引理 14.1：\(u,v\in\mathbb R^n\)，\(u^Tv\ge0\)，则 \(\|UVe\|\le2^{-3/2}\|u+v\|^2\)。

证明思路：对 \(\alpha\beta\ge0\) 有 \(\sqrt{|\alpha\beta|}\le\frac12|\alpha+\beta|\)（算术–几何平均）。把指标分成 \(u_iv_i\ge0\) 与 \(<0\) 两组，由 \(u^Tv\ge0\)，前一组的和不小于后一组，于是 \(\|UVe\|\le\sqrt2\,\|[u_iv_i]_{\mathcal P}\|_1\le\sqrt2\sum_{\mathcal P}\frac14(u_i+v_i)^2\le2^{-3/2}\|u+v\|^2\)。∎

引理 14.2：若 \((x,\lambda,s)\in\mathcal N_{-\infty}(\gamma)\)，则 (14.11) 的步满足

证明思路：由 (14.11) 前两行，\(\Delta x\in\mathrm{Null}(A)\)，\(\Delta s=-A^T\Delta\lambda\)，故 \(\Delta x^T\Delta s=0\)（14.35）。第三行乘 \((XS)^{-1/2}\)，记 \(D=X^{1/2}S^{-1/2}\)：\(D^{-1}\Delta x+D\Delta s=(XS)^{-1/2}(-XSe+\sigma\mu e)\)。对 \(u=D^{-1}\Delta x\)、\(v=D\Delta s\)（\(u^Tv=0\)）用引理 14.1，再利用 \(x_is_i\ge\gamma\mu\) 估计 \(\sum1/(x_is_i)\le n/(\gamma\mu)\)，得到界 \(2^{-3/2}(1-2\sigma+\sigma^2/\gamma)n\mu\le2^{-3/2}(1+1/\gamma)n\mu\)。∎

定理 14.3：存在与 \(n\) 无关的常数 \(\delta\) 使

证明思路。第一步，证明步长下界 \(\alpha_k\ge2^{3/2}\frac{\sigma_k}n\gamma\frac{1-\gamma}{1+\gamma}\)：展开 \(x_i(\alpha)s_i(\alpha)=x_is_i+\alpha(x_i\Delta s_i+s_i\Delta x_i)+\alpha^2\Delta x_i\Delta s_i\)，用第三行方程和引理 14.2 给出下界，与邻域条件 \(x_i(\alpha)s_i(\alpha)\ge\gamma\mu(\alpha)\) 比较即得。第二步，利用 \(\Delta x^T\Delta s=0\)，

\(\sigma(1-\sigma)\) 在区间端点取最小值，取 \(\delta=2^{3/2}\gamma\frac{1-\gamma}{1+\gamma}\min\{\sigma_{\min}(1-\sigma_{\min}),\sigma_{\max}(1-\sigma_{\max})\}\) 即可。∎

推导拆解：(14.43) 中 \(\mu_{k+1}=(1-\alpha(1-\sigma))\mu_k\) 这个等式值得自己推一遍，它也解释了 14.7.2 日志里 \(\mu\) 的下降速度。 第一步，展开新点的内积：\((x+\alpha\Delta x)^T(s+\alpha\Delta s)=x^Ts+\alpha(s^T\Delta x+x^T\Delta s)+\alpha^2\Delta x^T\Delta s\)。 第二步，把 (14.11) 第三行 \(s_i\Delta x_i+x_i\Delta s_i=-x_is_i+\sigma\mu\) 对 \(i\) 求和：\(s^T\Delta x+x^T\Delta s=-x^Ts+n\sigma\mu=-(1-\sigma)n\mu\)。 第三步，\(\Delta x^T\Delta s=0\)（\(\Delta x\) 在 \(A\) 的零空间里，\(\Delta s=-A^T\Delta\lambda\)，所以 \(\Delta x^T\Delta s=-(A\Delta x)^T\Delta\lambda=0\)）。 合起来 \(n\mu_{k+1}=n\mu_k-\alpha(1-\sigma)n\mu_k\)。读法：走全步 \(\alpha=1\) 时，间隙恰好按 \(\sigma\) 的比例缩小；\(\sigma=0.01\) 就是每步缩小 100 倍。

定理 14.4（复杂度）：若初始点在 \(\mathcal N_{-\infty}(\gamma)\) 中且 \(\mu_0\le1/\epsilon^\kappa\)，则存在 \(K=O(n\log1/\epsilon)\) 使 \(k\ge K\) 时 \(\mu_k\le\epsilon\)。

证明：对 (14.37) 取对数递推，用 \(\log(1+\beta)\le\beta\)：\(\log\mu_k\le-k\delta/n+\kappa\log\frac1\epsilon\)。要使右端 \(\le\log\epsilon\)，只需 \(k\ge K=\frac{(1+\kappa)n}\delta\log\frac1\epsilon\)。∎（原书证明中此处误印为 \(n\omega\)，应为 \(n\)。）

这个界是理论保证；实际算法（下一节）的迭代次数远少于 \(O(n)\)，基本与规模无关。

金融直觉：(14.37) 说 \(\mu\) 每步至少按固定比例 \(\delta/n\) 下降，这是几何（复利式）衰减，和按固定比例摊销的贷款余额一样。\(k\) 步后 \(\mu_k\le(1-\delta/n)^k\mu_0\approx e^{-k\delta/n}\mu_0\)，要降到 \(\epsilon\) 只需 \(k\approx\frac n\delta\ln\frac{\mu_0}\epsilon\)。关键在 \(\ln\frac1\epsilon\)：精度每提高 10 倍，只多花一个固定的步数，正如资金翻倍所需的年数（72 法则）不依赖本金大小。单纯形法的最坏情况则是"顶点个数"那样的指数增长。

14.3.2 其他原始–对偶方法（简介）

• 短步路径跟踪：取略小于 1 的保守 \(\sigma\)，使全步 \(\alpha=1\) 不离开严格的 \(\mathcal N_2\) 邻域。理论漂亮（\(O(\sqrt n\log1/\epsilon)\)），但进展慢。
• Mizuno–Todd–Ye 预测–校正法：使用两个嵌套的 \(\mathcal N_2\) 邻域，预测步（\(\sigma=0\)）从内邻域走到外邻域边界以大幅降低 \(\mu\)，校正步（\(\sigma=1\)，\(\alpha=1\)）拉回内邻域。\(\mu_k\) 超线性收敛到零。
• 势函数下降法：不显式跟踪中心路径，而用 Tanabe–Todd–Ye 势函数

要求每步固定量下降。第一项推动对偶间隙下降，第二项阻止单个乘积独立趋零。取 \(\sigma_k\equiv n/(n+\sqrt n)\) 即可保证 \(\Phi_\rho\) 每步固定下降。

14.4 实用算法：Mehrotra 预测–校正

几乎所有通用 LP 内点代码都基于 Mehrotra 预测–校正算法（predictor–corrector），它有两个特点：

1. 在方向上加一个校正步（corrector），补偿线性化忽略的二阶项，更紧密地跟踪通往解的轨迹；
2. 自适应地选择中心化参数 \(\sigma\)：先试算纯牛顿方向，若它能大幅降低 \(\mu\)，就少中心化；否则多中心化。

具体步骤（每步只做一次矩阵分解，三个右端共用）：

1. 预测步：在 (14.15) 中取 \(\sigma=0\)，

2. 沿预测方向不违反非负性的最大步长：

（原书 (14.22) 中第二个因子的步长也印成 \(\alpha^{\rm pri}_{\rm aff}\)，按 Mehrotra 原意应为对偶步长。）

3. 自适应中心化：\(\sigma=(\mu_{\rm aff}/\mu)^3\)。预测方向效果好时 \(\mu_{\rm aff}\ll\mu\)，\(\sigma\) 接近 0。

4. 校正 + 中心化：真实的互补条件是 \((x_i+\Delta x_i)(s_i+\Delta s_i)=\sigma\mu\)，线性化丢掉了二阶项 \(\Delta x_i\Delta s_i\)。用预测步的 \(\Delta x^{\rm aff}_i\Delta s^{\rm aff}_i\) 近似补上：

5. 原始与对偶分别取步长：

\(\eta\in[0.9,1)\)，接近解时 \(\eta\to1\) 以加速收敛。\(x\) 用原始步长更新，\((\lambda,s)\) 用对偶步长更新。

金融直觉：Mehrotra 的"预测 + 校正"和"久期 + 凸性"的思路一样。预测步是一阶近似，相当于只用久期估算债券价格变化；它丢掉了二阶项 \(\Delta x_i\Delta s_i\)。校正步用预测步算出来的 \(\Delta x^{\rm aff}_i\Delta s^{\rm aff}_i\) 当作二阶项的估计补回去，相当于加上凸性修正。两者用的是同一个系数矩阵，只换右端，所以只需一次矩阵分解，额外成本很低。 \(\sigma=(\mu_{\rm aff}/\mu)^3\) 的直觉是"看路况决定车速"：试走一步发现间隙能降得很多（\(\mu_{\rm aff}/\mu\) 小），说明前路通畅，就少做中心化；降得不多，说明离边界太近，就多做中心化。三次方是经验选择，没有理论推导。 原始、对偶分别取步长，是因为 \(x\) 与 \(s\) 离各自边界的距离不同，强制同一步长会让离边界远的一方白白少走。

注意：上述形式的 Mehrotra 算法没有收敛理论，甚至能构造出发散的例子。可以加保护措施纳入 14.3 节的收敛框架，但多数软件不这么做，因为它在实践中表现极好。现代软件还常加 Gondzio 的"高阶校正"。

14.5 每步的线性代数

主要计算量在于解 (14.15)、(14.20)、(14.23)。系数矩阵大而稀疏，但可以化为更紧凑的对称形式。

增广系统（augmented system）：由于 \(x,s>0\)，从第三行消去 \(\Delta s\)：

正规方程（normal equations）：再消去 \(\Delta x\)，

之所以叫"正规方程"，是因为 (14.28a) 可看作系数矩阵为 \(DA^T\) 的最小二乘问题的正规方程（本册第 10 章 10.2.3 节）。多数实现基于 (14.28)：对 \(AD^2A^T\) 做稀疏 Cholesky 分解，再回代恢复 \(\Delta s,\Delta x\)。

推导拆解：从 (14.15) 消元得到 (14.28)，记第三行右端为 \(-r_{xs}\)，\(r_{xs}=XSe-\sigma\mu e\)。 第一步，由第一行解出 \(\Delta s=-r_c-A^T\Delta\lambda\)，这就是 (14.28b)。 第二步，第三行 \(S\Delta x+X\Delta s=-r_{xs}\)，两边左乘 \(S^{-1}\)（\(s>0\) 所以可逆）：\(\Delta x=-S^{-1}r_{xs}-S^{-1}X\Delta s\)。注意 \(S^{-1}r_{xs}=x-\sigma\mu S^{-1}e\)，这就是 (14.28c)。 第三步，把第一步的 \(\Delta s\) 代入第二步，得 \(\Delta x=-S^{-1}r_{xs}+D^2r_c+D^2A^T\Delta\lambda\)（\(D^2=S^{-1}X\)）。 第四步，代入第二行 \(A\Delta x=-r_b\)，把含 \(\Delta\lambda\) 的项留在左边，就得到 (14.28a)。 和读者熟悉的回归对照：\(AD^2A^T\) 的形状与加权最小二乘（WLS）的 \(X^TWX\) 相同，权重 \(D^2_{ii}=x_i/s_i\)。接近最优时，"该为正"的变量 \(x_i/s_i\to\infty\)，"该为零"的变量 \(x_i/s_i\to0\)，权重悬殊越来越大，这正是下面说的后期病态。 Cholesky 分解是把对称正定矩阵写成 \(LL^T\)（\(L\) 下三角）的方法，之后解方程只需两次回代，可以理解为对称矩阵专用的、更快更稳的"求逆"。

两个需要注意的数值问题：

• 后期病态。接近解时，\(x_i\to0\)、\(s_i>0\) 的分量使 \(D^2_{ii}\to0\)，反之 \(D^2_{ii}\to\infty\)，\(AD^2A^T\) 极度病态甚至数值奇异。通用 Cholesky 需要小改动（如遇到极小主元时置为大数）。奇妙的是，由于结构特殊，算出的方向仍然足够好。
• 稠密列。\(A\) 中只要有一列稠密，\(AD^2A^T\) 就会完全稠密。增广系统 (14.26) 没有这个问题，也能自然处理自由变量，但需要稀疏对称不定分解，更复杂。

量化提示：组合优化里的预算约束 \(\mathbf 1^Tw=1\)、因子暴露约束 \(B^Tw=0\) 都是稠密行；情景型模型（CVaR、MAD、\(\ell_1\) 跟踪）的情景收益矩阵则整块稠密。这正是 14.7 节实验中内点法每步代价偏高的原因。自己写求解器时，要么用增广系统，要么利用 Sherman–Morrison–Woodbury 公式（附录 A）把稠密的低秩部分单独处理。

14.6 推广（简介）

原始–对偶框架可以推广到：

• 单调线性互补问题（LCP）：求 \(x,s\) 使 \(s=Mx+q\)，\((x,s)\ge0\)，\(x^Ts=0\)，\(M\) 半正定。美式期权定价的离散化（障碍问题）就是 LCP。
• 凸二次规划 \(\min c^Tx+\frac12x^TGx\) s.t. \(Ax=b,\ x\ge0\)，\(G\succeq0\)：第 16b 章 16.9 节（原书 §16.7）详细讨论，均值–方差组合优化即属此类。
• 非线性规划：把 KKT 写成类似 (14.3) 的形式，用牛顿法并截短步长保持正性，第 17 章 17.3 节（原书 §17.2）讨论。
• 半定规划（SDP）：变量为须半正定的对称矩阵。最近相关阵修复、某些稳健组合问题是 SDP。

LCP 与凸 QP 的推广保留了 LP 算法的收敛性与多项式复杂度。

14.7 量化实战

14.7.1 在哪里用

• 读懂求解器日志。MOSEK、Gurobi barrier、CPLEX barrier、HiGHS-IPM 的每行日志通常打印原始残差 \(\|r_b\|\)、对偶残差 \(\|r_c\|\)、对偶间隙或 \(\mu\)。理解本章后你会知道：残差先降到零（全步之后保持），\(\mu\) 随后以几乎固定的比例下降；若 \(\mu\) 停滞而残差不降，多半是问题不可行或数值病态。
• 对偶间隙做停止准则。内点法天然给出对偶界，可以量化"当前组合离最优还有多远"，适合盘中需要时效的再优化：间隙小于预期交易成本即可提前停止。
• 顶点解与交叉（crossover）。内点法给出的是内部近似解，若需要精确的顶点（稀疏）解，求解器会做 crossover。对需要严格稀疏持仓或要拿基状态做后续分析的场景，要打开这一选项或改用单纯形。

14.7.2 代码一：自己实现 Mehrotra 预测–校正

按原书习题 14.12 的方法构造已知解的随机 LP：\(x^*\) 前 \(m\) 个分量为正、其余为零，\(s^*\) 相反，\(c=A^T\lambda^*+s^*\)，\(b=Ax^*\)。

import numpy as np
from scipy.optimize import linprog

def mehrotra_lp(A, b, c, tol=1e-9, eta=0.99, max_iter=100, verbose=True):
    """Algorithm 14.3：Mehrotra 预测-校正原始-对偶内点法，标准形 LP。
    线性系统用正规方程 (14.28)：A D^2 A^T dlam = rhs，D^2 = X S^{-1}。"""
    m, n = A.shape
    x, lam, s = np.ones(n) * 10, np.zeros(m), np.ones(n) * 10   # 只要求 x,s>0
    def step_len(v, dv):
        neg = dv < 0
        return min(1.0, np.min(-v[neg] / dv[neg])) if neg.any() else 1.0
    def solve(rc, rb, rxs):
        # 解 [0 A^T I; A 0 0; S 0 X][dx;dl;ds] = [-rc; -rb; -rxs]
        d2 = x / s
        M = (A * d2) @ A.T
        rhs = -rb + A @ (-d2 * rc + rxs / s)
        dl = np.linalg.solve(M, rhs)
        ds = -rc - A.T @ dl
        dx = -(rxs + x * ds) / s
        return dx, dl, ds
    for k in range(max_iter):
        rb, rc = A @ x - b, A.T @ lam + s - c
        mu = x @ s / n
        if verbose:
            print(f"{k:2d}  mu={mu:9.2e}  |rb|={np.linalg.norm(rb):8.1e}  "
                  f"|rc|={np.linalg.norm(rc):8.1e}  obj={c@x:12.6f}")
        if mu < tol and np.linalg.norm(rb) < tol * (1 + np.linalg.norm(b)) \
                and np.linalg.norm(rc) < tol * (1 + np.linalg.norm(c)):
            return x, lam, s, k
        # 预测步 (14.20)
        dxa, dla, dsa = solve(rc, rb, x * s)
        ap, ad = step_len(x, dxa), step_len(s, dsa)
        mu_aff = (x + ap * dxa) @ (s + ad * dsa) / n          # (14.22)
        sigma = (mu_aff / mu) ** 3                            # 自适应中心化
        # 校正 + 中心化 (14.23)
        dx, dl, ds = solve(rc, rb, x * s + dxa * dsa - sigma * mu)
        ap = min(1.0, eta * step_len(x, dx))                  # (14.25)
        ad = min(1.0, eta * step_len(s, ds))
        x, lam, s = x + ap * dx, lam + ad * dl, s + ad * ds   # 原始/对偶分别步长
    raise RuntimeError("not converged")

# 按习题 14.12 构造已知解的随机 LP
rng = np.random.default_rng(0)
m, n = 50, 120
A = rng.standard_normal((m, n))
xs = np.zeros(n); xs[:m] = rng.uniform(0.5, 2, m)
ss = np.zeros(n); ss[m:] = rng.uniform(0.5, 2, n - m)
lam_s = rng.standard_normal(m)
c, b = A.T @ lam_s + ss, A @ xs
x, lam, s, it = mehrotra_lp(A, b, c)
print("迭代次数:", it, " ||x - x*|| =", np.linalg.norm(x - xs),
      " 对偶间隙 c'x-b'lam =", c @ x - b @ lam)
r = linprog(c, A_eq=A, b_eq=b, bounds=[(0, None)] * n, method="highs-ipm")
print("HiGHS-IPM 目标:", r.fun, " 自写目标:", c @ x, " HiGHS 迭代:", r.nit)

输出：

 0  mu= 1.00e+02  |rb|= 7.9e+02  |rc|= 1.2e+02  obj= 1373.510406
 1  mu= 2.42e+01  |rb|= 9.6e+01  |rc|= 3.3e+01  obj=  632.328463
 2  mu= 3.39e+00  |rb|= 9.6e-01  |rc|= 3.3e-01  obj=  344.738557
 3  mu= 8.47e-01  |rb|= 2.2e-01  |rc|= 8.3e-02  obj=   97.131765
 4  mu= 4.19e-01  |rb|= 1.0e-01  |rc|= 2.3e-02  obj=   66.559230
 5  mu= 2.43e-01  |rb|= 1.3e-02  |rc|= 6.7e-03  obj=   51.791762
 6  mu= 5.45e-02  |rb|= 2.0e-03  |rc|= 1.5e-03  obj=   37.154846
 7  mu= 1.32e-02  |rb|= 4.0e-04  |rc|= 3.1e-04  obj=   34.285874
 8  mu= 9.04e-04  |rb|= 3.3e-05  |rc|= 1.1e-05  obj=   33.516065
 9  mu= 9.53e-06  |rb|= 3.4e-07  |rc|= 1.2e-07  obj=   33.441618
10  mu= 9.53e-08  |rb|= 3.4e-09  |rc|= 1.2e-09  obj=   33.440846
11  mu= 9.53e-10  |rb|= 3.4e-11  |rc|= 1.2e-11  obj=   33.440838
迭代次数: 11  ||x - x*|| = 9.56524778620355e-08  对偶间隙 c'x-b'lam = 1.1435405156134948e-07
HiGHS-IPM 目标: 33.44083822507707  自写目标: 33.44083830303802  HiGHS 迭代: 9

这就是一份"求解器日志"。前几步原始、对偶残差和 \(\mu\) 同步下降；第 8 步以后每步 \(\mu\) 缩小约 100 倍——步长 \(\eta=0.99\) 接近全步、\(\sigma\) 接近 0 时，\(\mu\) 每步乘以约 \(1-\eta\)。从远离可行的起点出发，11 次迭代就把对偶间隙降到 \(10^{-7}\)，与 HiGHS 的结果一致。

14.7.3 代码二：\(\ell_1\) 指数跟踪——单纯形法与内点法的迭代次数

指数增强 / 被动跟踪中，常用 MAD（平均绝对偏差）度量跟踪误差：

把每期误差 \(e_t\) 拆成 \(u_t-v_t\)（第 13 章技巧），就是一个 LP。

import numpy as np, time
from scipy.optimize import linprog
from scipy.sparse import hstack, vstack, eye, csr_matrix

def l1_tracking(n, T, seed=1):
    """指数跟踪：min (1/T) sum_t |R_t w - r_idx,t|, s.t. sum w = 1, 0<=w<=0.05
    拆分 e_t = u_t - v_t, u,v>=0（第 13 章的 x = x+ - x- 技巧）"""
    rng = np.random.default_rng(seed)
    f = 0.01 * rng.standard_normal(T)
    R = np.outer(f, rng.uniform(0.6, 1.4, n)) + 0.02 * rng.standard_normal((T, n))
    idx_w = rng.dirichlet(np.full(n, 0.2))   # 头部集中的指数权重
    r_idx = R @ idx_w
    c = np.r_[np.zeros(n), np.full(2 * T, 1.0 / T)]
    A_eq = vstack([hstack([csr_matrix(R), -eye(T), eye(T)]),
                   csr_matrix(np.r_[np.ones(n), np.zeros(2 * T)][None, :])])
    b_eq = np.r_[r_idx, 1.0]
    bounds = [(0, 0.05)] * n + [(0, None)] * (2 * T)
    out = {}
    for meth in ["highs-ds", "highs-ipm"]:
        t0 = time.perf_counter()
        r = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method=meth)
        out[meth] = (r.nit, time.perf_counter() - t0, r.fun)
    return out

print(f"{'n':>5}{'T':>6} | {'单纯形迭代':>8}{'用时s':>8} | {'内点迭代':>8}{'用时s':>8} | 跟踪误差(日均|e|)")
for n, T in [(50, 250), (200, 500), (500, 1000)]:
    o = l1_tracking(n, T)
    print(f"{n:5d}{T:6d} | {o['highs-ds'][0]:10d}{o['highs-ds'][1]:8.2f} | "
          f"{o['highs-ipm'][0]:10d}{o['highs-ipm'][1]:8.2f} | {o['highs-ipm'][2]*1e4:.2f}bp")

输出：

    n     T |    单纯形迭代     用时s |     内点迭代     用时s | 跟踪误差(日均|e|)
   50   250 |        400    0.01 |         21    0.02 | 22.31bp
  200   500 |       1007    0.17 |         21    0.20 | 0.52bp
  500  1000 |       2999    2.26 |         23    1.15 | 1.96bp

（用时一栏随机器而异，迭代次数是确定的。）规律很清楚：单纯形迭代次数随规模大致线性增长（几百到几千），内点法稳定在 21–23 次。用时则取决于每步代价：本例情景收益矩阵 \(R\) 是稠密的，内点法每步的 \(AD^2A^T\) 分解很贵，所以小问题上单纯形更快、大问题上内点法反超。这也说明为什么商业求解器默认会"并发"运行单纯形和内点法，谁先完成用谁。跟踪误差一栏说明了个股 5% 上限的作用：头部权重超过 5% 的成分股无法完全复制，误差不为零。

本章小结

原始–对偶内点法把 LP 的 KKT 条件写成 \(F(x,\lambda,s)=0\) 加 \((x,s)\ge0\)，对等式部分用牛顿法，同时严格保持 \((x,s)>0\) 以避开伪解。关键概念是中心路径（所有 \(x_is_i\) 相等的严格可行曲线）、对偶度量 \(\mu=x^Ts/n\) 和中心化参数 \(\sigma\)：方向在纯牛顿（\(\sigma=0\)）和中心化（\(\sigma=1\)）之间权衡。路径跟踪方法让迭代点留在中心路径的锥形邻域里，长步算法每步 \(\mu\) 至少缩小 \(1-\delta/n\)，因此是多项式的。实用的 Mehrotra 预测–校正算法用预测步自适应选 \(\sigma=(\mu_{\rm aff}/\mu)^3\)，再加二阶校正，一次分解、原始对偶分别取步长，虽无收敛理论但实践中极好，迭代次数几乎不随规模增长。每步的核心是解正规方程 \(AD^2A^T\Delta\lambda=\cdots\)，要注意后期病态和稠密列。这一框架直接推广到凸二次规划（第 16 章）和非线性规划（第 17 章）。

练习

基础

1. 验证 14.2.1 节的伪解例子：写出 \(F\) 的三个分量并逐一代入。再找出该问题 \(F=0\) 的全部解。
2. 证明：(i) \(\theta_1<\theta_2\) 时 \(\mathcal N_2(\theta_1)\subset\mathcal N_2(\theta_2)\)；\(\gamma_2\le\gamma_1\) 时 \(\mathcal N_{-\infty}(\gamma_1)\subset\mathcal N_{-\infty}(\gamma_2)\)；(ii) \(\gamma\le1-\theta\) 时 \(\mathcal N_2(\theta)\subset\mathcal N_{-\infty}(\gamma)\)；(iii) \(\mathcal N_{-\infty}(1)=\mathcal N_2(0)=\mathcal C\)。 提示：(ii) \(|x_is_i-\mu|\le\|XSe-\mu e\|_2\le\theta\mu\)。
3. 证明 (14.11) 的系数矩阵非奇异当且仅当 \(A\) 行满秩（设 \(x,s>0\)）。
4. 证明 \(AD^2A^T\) 对称正定当且仅当 \(A\) 行满秩（\(D\) 对角元全正）。若 \(D\) 恰有 \(m\) 个正对角元、其余为零，结论是否仍成立？这和内点法后期病态有什么联系？
5. 在 14.7.2 的代码中把 \(\sigma\) 固定为 0（纯仿射尺度）和固定为 0.5，比较迭代次数。

进阶

6. 证明 Tanabe–Todd–Ye 势函数满足：若某个 \(x_is_i\to0\) 而 \(\mu\) 不趋于零，则 \(\Phi_\rho\to\infty\)。
7. 对含自由变量 \(y\) 的 LP \(\min c^Tx+d^Ty\) s.t. \(A_1x+A_2y=b\)，\(x\ge0\)，写出最优性条件、原始–对偶步方程和增广系统形式，解释为什么它不能化为对称正定的正规方程。
8. 对轨迹 \(\mathcal H\)：\(F(\hat x(\tau),\hat\lambda(\tau),\hat s(\tau))=(1-\tau)(r_c,r_b,XSe)\)，求 \(\tau=0\) 处一阶、二阶导数满足的方程，说明 Mehrotra 校正项 \(-\Delta X^{\rm aff}\Delta S^{\rm aff}e\) 正是二阶 Taylor 项。
9. 在 14.7.3 的跟踪问题中，分别固定 \(T\) 增大 \(n\)、固定 \(n\) 增大 \(T\)，记录内点法的迭代次数和每次迭代的平均用时。用正规方程 \(AD^2A^T\) 的维数（等式约束个数 \(T+1\)）和稠密程度解释你看到的规律。 提示：迭代次数几乎不变；每步代价主要由 \((T+1)\times(T+1)\) 稠密矩阵的分解决定。
10. 用 14.7.2 的代码求解一个不可行的 LP（如 \(x_1+x_2=-1,\ x\ge0\)），观察 \(\mu\) 与残差的行为，总结如何从日志判断不可行。

原书推荐习题：14.1（伪解）、14.2、14.5（中心路径与邻域）、14.9、14.11（正规方程与增广系统的适用条件）、14.10（预测–校正的 Taylor 展开）、14.12（亲手实现 Mehrotra 算法，最有实践价值）。原书第 14 章共 14.1–14.12 题。

原书对照

更正说明：原书 (14.22) 第二个因子的步长应为对偶步长；定理 14.4 证明中的 \(n\omega\) 应为 \(n\)。原书正文页码 = PDF 页码减 19（已用 PDF 页眉核对，第 12–16 章均如此）。