第 29 章 线性规划
本章对应原书第 29 章。线性规划(linear programming,LP)研究的是:在一组线性等式和不等式约束下,最大化或最小化一个线性函数。它是运筹学的基石,也是组合构建、指数复制、执行调度中最常用的建模工具之一。原书用代数方法讲单纯形算法(simplex algorithm):把问题写成"松弛型",每次做一次"转轴",像高斯消元一样不断改写方程,直到最优解一眼可见;然后用对偶性证明这个解确实最优,并顺手得到每条约束的"影子价格"。本章把这条主线讲透:先建模,再单纯形,再对偶,再处理初始可行解。第 04 册第 13 章从数值最优化的角度讲单纯形法的矩阵实现(修正单纯形法、LU 更新)和第 14 章的内点法,本章与之互补,侧重组合结构与证明。量化实战部分用 LP 做部分复制的指数跟踪、增强指数的 alpha-跟踪误差前沿,以及用 Farkas 引理检验无套利。
学习目标
读完本章,你应当能够:
- 把实际问题(包括最短路径、最大流、最小费用流、组合构建)写成线性规划,并化为标准型 \(\max c^Tx,\ Ax\le b,\ x\ge0\) 和松弛型。
- 手工执行单纯形算法:选入基变量、做最小比值检验、转轴、读出基本解;理解退化、循环和 Bland 规则。
- 写出任意线性规划的对偶,证明弱对偶,理解强对偶的构造性证明,并从单纯形的最终松弛型中读出对偶最优解(影子价格)。
- 用辅助线性规划(两阶段法)判定可行性、获得初始基本可行解,说出线性规划的三种结局。
- 用互补松弛条件检验最优性;用 Farkas 引理理解"无套利 ⟺ 存在正的状态价格"。
- 用
scipy.optimize.linprog求解指数跟踪 LP,解读约束的影子价格,并说明 L1 跟踪、基数约束、整数手数等建模选择的后果。
读前导读
这一章在解决什么问题
一句话:在一堆"不能超过""至少要有"的线性限制下,把一个线性目标做到最大或最小。你在 CPA 管理会计和 CFA 公司金融里其实已经做过线性规划,只是没有叫这个名字。
最贴近的例子是资本限额下的资本预算(capital rationing)。假设今年能投的资本只有 1 亿元,候选项目有 5 个,每个项目可以投一部分(比如按比例参与),每个项目有 NPV。教科书的做法是按"盈利指数 PI = NPV / 投资额"从高到低排序、一直投到钱用完。这正是只有一条约束的线性规划的最优解。但现实中约束往往不止一条:今年资本上限、明年资本上限、某个部门的人力上限、某类风险敞口上限……一旦约束有两条以上,"按 PI 排序"就失灵了,因为一个项目可能在资本上很划算、在人力上很昂贵。线性规划就是这类问题的通用解法:把"每个项目投多少"设成变量 \(x_j\),目标是 \(\max\sum_j\text{NPV}_jx_j\),每条资源限制写成一条不等式。
本章还有第二个、对金融人更重要的收获:对偶与影子价格。求出最优方案的同时,单纯形法会告诉你"每条约束放宽一单位,目标能多挣多少"。在资本预算里,这就是资本的影子成本——如果今年再多 100 万元预算能多带来 18 万元 NPV,那么任何融资成本低于 18% 的额外融资都值得争取。在组合构建里,它告诉你行业中性、单票上限这些约束各自"吃掉"了多少 alpha。最后,对偶理论的一个变形(Farkas 引理)正是"无套利 ⟺ 存在状态价格",也就是你在 CFA 衍生品里用单步二叉树定价时背后的那条定理。
与第 04 册的分工:第 04 册第 13 章 线性规划:单纯形法 用矩阵语言讲单纯形法怎样在计算机上高效、稳定地实现,第 14 章 内点法 讲另一类算法;第 12 章 约束优化理论 的 KKT 条件是本章互补松弛在非线性问题里的推广。本章用"手算方程组"的方式讲同一个算法,更容易看清每一步在干什么,适合先读。
需要先想起来的数学
1. 求和记号与矩阵记号。 \(\sum_{j=1}^n a_{ij}x_j\) 就是"第 \(i\) 行的系数和变量逐个相乘再相加"。把所有行叠起来就写成 \(Ax\),其中 \(A\) 是 \(m\) 行 \(n\) 列的系数表。\(c^Tx\)(\(T\) 表示转置,把列向量横过来)就是 \(\sum_jc_jx_j\),例如 \(c=(3,1,2)\)、\(x=(8,4,0)\) 时 \(c^Tx=24+4+0=28\)。\(A^Ty\) 是把系数表"转个方向"来乘,第 \(j\) 个分量是 \(\sum_ia_{ij}y_i\),即按列求和。向量不等式 \(x\ge0\) 表示每个分量都 \(\ge0\)。见 第 00 册第 06 章 线性代数速成 和 第 07 章 概率中的分析工具 的求和记号部分。
2. 解线性方程组(代入消元)。 单纯形法的每一步本质上就是中学的代入法:从一个方程里解出某个变量,代入其余方程。例如由 \(x_6=36-4x_1-x_2-2x_3\) 解出 \(x_1=9-\frac{x_2}{4}-\frac{x_3}{2}-\frac{x_6}{4}\),再代进别的式子。能熟练做这个,本章的计算就没有障碍。见 第 00 册第 06 章 的线性方程组部分。
3. 拉格朗日乘子与"影子价格"。 在带约束的最优化里,拉格朗日乘子 \(\lambda\) 的含义是"约束右端放宽一单位,最优值变化多少"。本章的对偶变量 \(y_i\) 就是线性规划里的拉格朗日乘子。你不需要会推拉格朗日条件,只要记住这个经济含义。见 第 00 册第 05 章 多元微积分与优化。
4. 读证明的几个词。 "当且仅当(⟺)"表示两个方向都成立;"恰有一个成立"表示两件事不会同时成立、也不会同时不成立;"引理"是为主定理服务的小结论;证明末尾的 \(\square\) 表示证毕。\(\forall\) 读作"对所有",\(\exists\) 读作"存在"。\(\arg\min\) 表示"让某个量最小的那个下标",而不是最小值本身。见 第 00 册第 08 章 读懂数学证明与符号。
怎么读这一章
核心必读是 29.1(建立几何直观)、29.2 的松弛型记号、29.4.2 的完整例题、29.5 的对偶与影子价格,以及 29.8.2–29.8.3 两个量化实战。建议的顺序是:先读 29.1 和 29.4.2,拿纸笔把三次转轴跟着算一遍,这比读任何形式化定义都有用;再回头读 29.2 的记号;然后读 29.5,重点是"怎样想到对偶"、弱对偶的两行证明和影子价格。
第一次可以只看结论的部分:29.3 中的最大流、多商品流(需要第 26 章的背景);29.4.3 的 PIVOT 伪代码和引理 29.1(它只是把手算过程写成通用公式);29.4.5 退化与循环的细节和引理 29.4–29.7;29.6 的两阶段法(只需记住"先求一个可行点,再求最优")。强对偶的证明(29.5.5)有一定难度,第一遍可以只接受结论"单纯形法顺手给出了对偶最优解"。
29.1 什么是线性规划
29.1.1 一个政治问题
原书的引例:一位候选人希望在城市、郊区、农村三类选区(登记选民分别为 100,000、200,000、50,000)中各赢得至少一半选票,即 50,000、100,000、25,000 张。可投放广告的议题有四个:修路、枪支管制、农业补贴、专用于公共交通的汽油税。每在某议题上花 1000 美元广告,在各类选区赢得(负数为失去)的选票数(千张)如下(原书图 29.1):
| 政策 | 城市 | 郊区 | 农村 |
|---|---|---|---|
| 修路 | −2 | 5 | 3 |
| 枪支管制 | 8 | 2 | −5 |
| 农业补贴 | 0 | 0 | 10 |
| 汽油税 | 10 | 0 | −2 |
试错方案"修路 20、枪支 0、补贴 4、汽油税 9(千美元)"在城市赢得 \(20(-2)+0+4(0)+9(10)=50\),郊区 \(20(5)=100\),农村 \(20(3)+0+4(10)+9(-2)=82\),花费 33 千美元。能更省吗?设 \(x_1,\dots,x_4\) 为四项广告支出(千美元):
(负面广告可能存在,但广告费不能是负的,所以要求非负。)29.8 节的代码求出最优解是 \(x=(2050/111,\ 425/111,\ 0,\ 625/111)\),最少花费 \(3100/111\approx27.93\) 千美元,比试错方案省约 15%。
29.1.2 一般定义
线性函数 \(f(x_1,\dots,x_n)=\sum_ja_jx_j\)。\(f=b\) 是线性等式,\(f\le b\)、\(f\ge b\) 是线性不等式,统称线性约束(linear constraints)。注意不允许严格不等式(习题 29.5-4:\(\max x\) s.t. \(x<1\) 没有最优解)。线性规划问题是在有限个线性约束下最小化或最大化一个线性函数。
29.1.3 几何直观
二维例子:
满足所有约束的点称为可行解(feasible solution);它们组成平面上的一个凸多边形,即可行域(feasible region)。要最大化的函数称为目标函数(objective function),取值称为目标值。画出等值线 \(x_1+x_2=z\)(斜率为 \(-1\) 的平行直线族),把它往目标增大的方向平移,直到即将离开可行域为止。本例最优解为 \(x_1=2,x_2=6\),目标值 8。
关键观察:最优解总可以在可行域的某个顶点上取到。等值线最后与可行域的交集要么是一个顶点,要么是一条边(边上各点目标值相同,端点也是顶点)。在 \(n\) 维空间中,每个约束定义一个半空间,它们的交称为单纯形(simplex),目标函数的等值面是超平面,由于凸性,最优解仍在某个顶点上。
白话解释:为什么最优解一定在"角"上?目标是线性的,意味着沿任何方向走,目标值都按固定速度变化,不会先升后降。所以如果你站在可行域内部或者一条边的中间,总能沿某个方向继续走、让目标不变差,一直走到撞上边界、再撞上另一条边界,最终停在一个角上。资本预算里的对应现象是:在资本约束下按比例参与项目,最优方案通常是"几个项目全投、几个项目不投、最多一两个项目部分投",而不会是"每个项目都投一点"——那些"全投或不投"正是一条条约束取等,构成了顶点。
("凸"的意思是:可行域里任取两点,连线上的点也都可行。这里每条线性约束切出一个"半边空间",半边空间的交集一定是凸的。原书把这个凸区域称为"单纯形",这是 CLRS 的用法;在多数数学书里,"单纯形"专指三角形、四面体这类最简单的凸体,读其他书时注意区分。)
29.1.4 单纯形算法的思路与算法全景
单纯形算法从可行域的某个顶点出发,每次沿一条边移到目标值不更小的相邻顶点,直到所有相邻顶点都不更好。因为可行域是凸的、目标是线性的,局部最优就是全局最优(29.5 节用对偶性严格证明)。原书用代数而非几何描述这一过程:把 LP 写成松弛型,"从一个顶点走到相邻顶点"就是把一个基本变量换成非基本变量,称为转轴(pivot)。
算法全景:
- 单纯形法(Dantzig,1947):实践中通常很快,但对精心构造的输入需要指数时间(29.4.7 节的 Klee-Minty 例子)。
- 椭球法(ellipsoid algorithm,Khachian 1979):第一个多项式时间算法,实践中很慢。
- 内点法(interior-point methods,Karmarkar 首创):穿过可行域内部,中间解不一定是顶点;对大规模问题与单纯形法一样快,有时更快(第 04 册第 14 章)。
- 整数线性规划(integer linear program):再要求变量取整数,仅判定可行性就是 NP 难的(习题 34.5-3),没有已知多项式算法。
记号约定:变量写作 \(x=(x_1,\dots,x_n)\),某个具体取值写作 \(\bar x\)。
29.2 标准型与松弛型
29.2.1 标准型
标准型(standard form):给定实数 \(c_j\)、\(b_i\)、\(a_{ij}\),
紧凑形式 \(\max c^Tx\),s.t. \(Ax\le b\),\(x\ge0\),用三元组 \((A,b,c)\) 表示。\(x\ge0\) 称为非负约束。
术语:满足所有约束的 \(\bar x\) 是可行解,否则是不可行解;目标值最大的可行解是最优解。没有可行解的 LP 是不可行的(infeasible);有可行解但目标值无上界的 LP 是无界的(unbounded)。可行域无界时,最优值仍可能有限(习题 29.1-9)。
29.2.2 化为标准型的四条规则
| 不是标准型的原因 | 处理 |
|---|---|
| 目标是最小化 | 目标系数取负,变为最大化 |
| 变量 \(x_j\) 没有非负约束 | 用 \(x_j'-x_j''\) 替换,\(x_j',x_j''\ge0\) |
| 等式约束 \(f=b\) | 拆成 \(f\le b\) 与 \(f\ge b\) |
| \(\ge\) 约束 | 两边乘 \(-1\) 变为 \(\le\) |
两个 LP 等价是指:一个的每个目标值为 \(z\) 的可行解,在另一个中都有目标值为 \(z\)(一个最小化、一个最大化时为 \(-z\))的可行解,反之亦然。
例 \(\min\ -2x_1+3x_2\),s.t. \(x_1+x_2=7\),\(x_1-2x_2\le4\),\(x_1\ge0\)。取负得 \(\max\ 2x_1-3x_2\);\(x_2\) 无符号约束,替换为 \(x_2'-x_2''\);等式拆成两个不等式,\(\ge\) 的那个取负;把 \(x_2',x_2''\) 重命名为 \(x_2,x_3\),得标准型
29.2.3 松弛型
单纯形法希望除非负约束外都是等式。对不等式 \(\sum_ja_{ij}x_j\le b_i\),引入松弛变量(slack variable)
它度量第 \(i\) 个约束离"取等"还有多远。上例引入 \(x_4,x_5,x_6\),并用 \(z\) 表示目标值,略去"max""s.t."和显式的非负约束,得到松弛型(slack form):
等式左边的变量称为基本变量(basic variables),右边的称为非基本变量(nonbasic variables)。每个基本变量恰好出现在一个等式的左边;所有等式右边、以及目标函数中只出现非基本变量。
紧凑记号:\(N\) 为非基本变量的下标集,\(B\) 为基本变量的下标集,\(|N|=n\),\(|B|=m\)。松弛型用六元组 \((N,B,A,b,c,\nu)\) 表示:
注意减号:\(a_{ij}\) 是松弛型中"看上去"的系数的相反数。下标不必连续,由 \(B\)、\(N\) 决定。
白话解释:松弛变量就是"剩余额度"。资本预算里,约束"各项目投资合计 ≤ 1 亿"对应的松弛变量就是"没花掉的预算"。原约束 \(\le\) 换成"剩余额度 = 上限 − 已用额度,且剩余额度 ≥ 0",信息一点没丢,但不等式变成了等式,就可以像解方程组一样代入消元。
"基本变量 / 非基本变量"只是一种记账方式:把一部分变量放在等号左边(由别人决定),其余放在右边(自己可以自由设定)。令右边的非基本变量全取 0,左边就直接读出数值——这就是下一节的"基本解"。一开始,左边是全部松弛变量、右边是全部决策变量,对应"一个项目都不投、预算全部剩着"的起点。单纯形法每一步做的事,就是把某个决策变量从右边挪到左边、同时把一个松弛变量从左边挪到右边,即"开始投某个项目,直到某项预算被用完"。
例 松弛型
中 \(B=\{1,2,4\}\),\(N=\{3,5,6\}\),\(b=(8,4,18)\),\(c=(-1/6,-1/6,-2/3)\),\(\nu=28\),例如 \(a_{13}=-1/6\)(因为 \(x_1\) 的方程中是 \(+x_3/6\))。这个松弛型稍后会作为单纯形法例题的终点再次出现。
29.3 把问题写成线性规划
能识别"这个问题可以写成 LP"非常有用:一旦写成多项式规模的 LP,就可以用椭球法或内点法在多项式时间内求解,也可以直接交给现成求解器。
29.3.1 最短路径
单对最短路径(求 \(s\) 到 \(t\) 的 \(d_t\)):
为什么是最大化?若最小化,令所有 \(d_v=0\)(或更小)就满足约束,却没有求出最短路。最短路问题的解满足 \(d_v=\min_{(u,v)\in E}\{d_u+w(u,v)\}\)——这是满足约束的最大值。这与第 24 章的差分约束是同一个结构;约束 \(d_v-d_u\le w(u,v)\) 就是差分约束。
29.3.2 最大流
\(|V|^2\) 个变量;只为实际存在的边设变量时只需 \(O(V+E)\) 个约束(习题 29.2-5)。它的对偶就是最小割(习题 29.4-3)。
29.3.3 最小费用流
LP 求解已有专用算法的问题(Dijkstra、推送-重贴标签),通常不如专用算法快。LP 真正的威力在于新问题和没有已知高效算法的变体。最小费用流(minimum-cost flow)是最大流的推广:每条边除容量 \(c(u,v)\) 外还有费用 \(a(u,v)\),要从 \(s\) 送 \(d\) 单位流到 \(t\),使总费用最小:
例(原书图 29.3) 边(容量,费用):\(s\to x\)(5,2),\(s\to y\)(2,5),\(x\to y\)(1,3),\(x\to t\)(2,7),\(y\to t\)(4,1),送 4 单位。最优解 \(f_{sx}=2,f_{sy}=2,f_{xy}=1,f_{xt}=1,f_{yt}=3\),总费用 \(2\cdot2+5\cdot2+3\cdot1+7\cdot1+1\cdot3=27\)(已用 linprog 核对)。最小费用流可用于订单路由、资金调拨等"既有容量又有成本"的网络问题(第 26 章)。
29.3.4 多商品流
\(k\) 种商品共用一个网络,商品 \(i\) 用 \((s_i,t_i,d_i)\) 描述。各商品各自守恒,聚合流不超过边容量:
这是一个只问可行性的"空目标"LP。目前已知的唯一多项式时间算法就是把它写成 LP——LP 作为"通用求解器"的价值由此可见。多资产、多场所的订单拆分与路由在结构上就是多商品流。
29.4 单纯形算法
单纯形法可以看作"不等式版的高斯消元":高斯消元每步把方程组改写成结构更好的等价形式;单纯形法每步把松弛型改写成等价的松弛型,直到最优解一眼可见。
29.4.1 基本解
令所有非基本变量为 0,由等式算出基本变量的值,得到的解称为基本解(basic solution)。基本解中 \(\bar x_i=b_i\)(\(i\in B\))。若基本解可行(所有 \(b_i\ge0\)),称为基本可行解(basic feasible solution)。几何上,基本可行解就是可行域的顶点。
29.4.2 完整例题
松弛型:
基本解 \((0,0,0,30,24,36)\),目标值 0。对给定的非基本变量取值,若某约束的基本变量恰为 0,称该约束是紧的(tight)。
第 1 次转轴 \(x_1\) 在目标中系数为正,增大它能增大 \(z\)。增大到多少?\(x_4,x_5,x_6\) 分别在 \(x_1>30\)、\(12\)、\(9\) 时变负,第三个约束最紧,所以 \(x_1\) 最多增到 9,此时 \(x_6=0\)。交换 \(x_1\) 与 \(x_6\) 的角色:由 \(x_6\) 的方程解出 \(x_1=9-\frac{x_2}{4}-\frac{x_3}{2}-\frac{x_6}{4}\),代入其余方程:
\(x_1\) 称为入基变量(entering variable),\(x_6\) 称为出基变量(leaving variable)。新基本解 \((9,0,0,21,6,0)\),目标值 27。改写不改变 LP 本身:原基本解 \((0,0,0,30,24,36)\) 代入新方程仍成立,目标值 \(27-\frac34\cdot36=0\)。
推导拆解:这次转轴的代数可以逐步复原。
第一步,从 \(x_6\) 的方程 \(x_6=36-4x_1-x_2-2x_3\) 解出 \(x_1\):把 \(4x_1\) 移到左边、\(x_6\) 移到右边,得 \(4x_1=36-x_2-2x_3-x_6\),两边除以 4,得 \(x_1=9-\frac{x_2}{4}-\frac{x_3}{2}-\frac{x_6}{4}\)。
第二步,把它代入目标:\(z=3x_1+x_2+2x_3=3\big(9-\frac{x_2}{4}-\frac{x_3}{2}-\frac{x_6}{4}\big)+x_2+2x_3=27+\big(1-\frac34\big)x_2+\big(2-\frac32\big)x_3-\frac34x_6\),即 \(27+\frac{x_2}{4}+\frac{x_3}{2}-\frac{3x_6}{4}\)。
第三步,代入 \(x_4\) 的方程:\(x_4=30-x_1-x_2-3x_3=30-9+\frac{x_2}{4}+\frac{x_3}{2}+\frac{x_6}{4}-x_2-3x_3=21-\frac{3x_2}{4}-\frac{5x_3}{2}+\frac{x_6}{4}\)。\(x_5\) 同理:\(24-2\cdot9=6\),\(x_2\) 的系数 \(-2+\frac12=-\frac32\),\(x_3\) 的系数 \(-5+1=-4\),\(x_6\) 的系数 \(+\frac12\)。
读目标函数的方法:常数 27 是当前方案的目标值;\(+\frac{x_2}{4}\) 表示"再让 \(x_2\) 增加 1,目标还能多 \(\frac14\)";\(-\frac{3x_6}{4}\) 表示"如果把第三条约束留出 1 单位余量(\(x_6\) 从 0 变 1),目标会少 \(\frac34\)"。这些系数在运筹学里叫检验数或降低成本(reduced cost),作用相当于资本预算里"再多投一单位该项目的边际 NPV"。只要还有正的检验数,就还有改进空间。
第 2 次转轴 现在不能增大 \(x_6\)(系数为负,会降低目标)。选 \(x_3\):三个约束分别限制 \(x_3\le18\)、\(42/5\)、\(3/2\),第三个最紧,\(x_3\) 入基、\(x_5\) 出基:
基本解 \((33/4,0,3/2,69/4,0,0)\),目标值 \(111/4\)。
第 3 次转轴 只能增大 \(x_2\)。三个约束给出的上界为 132、4、\(\infty\)(\(x_4\) 的方程中 \(x_2\) 系数为正,增大 \(x_2\) 反而使 \(x_4\) 增大,不构成限制)。\(x_2\) 入基、\(x_3\) 出基:
目标函数中所有系数都为负——任何非基本变量从 0 增大都只会降低 \(z\)——所以当前基本解 \((8,4,0,18,0,0)\) 最优,目标值 28。回到原问题:\(x_1=8,x_2=4,x_3=0\),\(3\cdot8+4=28\)。
松弛变量的含义:\(x_4=18\) 表示第一个约束左边 \(8+4+0=12\) 比右边 30 少 18;\(x_5=x_6=0\) 表示后两个约束取等(是紧的)。
常见误区:即使输入全是整数,中间松弛型的系数和中间解也可能是分数;最终最优解也未必是整数——本例恰好是整数纯属巧合。这正是整数规划困难的根源之一。
29.4.3 转轴过程 PIVOT
输入松弛型 \((N,B,A,b,c,\nu)\)、出基下标 \(l\)、入基下标 \(e\),输出新松弛型:
def PIVOT(N, B, A, b, c, v, l, e):
# 1) 新基本变量 x_e 的方程:把 x_l 那一行解出 x_e
b_hat[e] = b[l] / A[l][e]
for j in N - {e}:
A_hat[e][j] = A[l][j] / A[l][e]
A_hat[e][l] = 1 / A[l][e]
# 2) 其余约束:把 x_e 的新表达式代入
for i in B - {l}:
b_hat[i] = b[i] - A[i][e] * b_hat[e]
for j in N - {e}:
A_hat[i][j] = A[i][j] - A[i][e] * A_hat[e][j]
A_hat[i][l] = -A[i][e] * A_hat[e][l]
# 3) 目标函数同样代入
v_hat = v + c[e] * b_hat[e]
for j in N - {e}:
c_hat[j] = c[j] - c[e] * A_hat[e][j]
c_hat[l] = -c[e] * A_hat[e][l]
# 4) 更新基/非基集合
N_hat = (N - {e}) | {l}
B_hat = (B - {l}) | {e}
return N_hat, B_hat, A_hat, b_hat, c_hat, v_hat
时间 \(O(mn)\)。只在 \(a_{le}>0\) 时调用,所以不会除以零。
引理 29.1 设 \(\bar x\) 为 PIVOT 之后的基本解,则 \(\bar x_j=0\)(\(j\in\hat N\)),\(\bar x_e=b_l/a_{le}\),\(\bar x_i=b_i-a_{ie}\hat b_e\)(\(i\in\hat B-\{e\}\))。
29.4.4 形式化的单纯形算法
def SIMPLEX(A, b, c):
N, B, A, b, c, v = INITIALIZE_SIMPLEX(A, b, c) # 29.6 节
while any(c[j] > 0 for j in N):
e = 选取某个 j in N 且 c[j] > 0 # 按某确定性规则
for i in B:
delta[i] = b[i] / A[i][e] if A[i][e] > 0 else INF
l = argmin_{i in B} delta[i] # 最小比值检验(ratio test)
if delta[l] == INF:
return "unbounded"
N, B, A, b, c, v = PIVOT(N, B, A, b, c, v, l, e)
return [b[i] if i in B else 0 for i in 1..n]
最小比值检验(minimum ratio test)选出对 \(x_e\) 增大限制最严的约束,其基本变量出基;若没有任何约束限制 \(x_e\)(所有 \(a_{ie}\le0\)),LP 无界。
引理 29.2 若 INITIALIZE-SIMPLEX 返回的松弛型基本解可行,则 SIMPLEX 返回的解可行;若返回"unbounded",LP 确实无界。
证明(三部分循环不变式:每次迭代开始时,(1) 松弛型与初始松弛型等价;(2) 对所有 \(i\in B\),\(b_i\ge0\);(3) 基本解可行)。关键是 (2) 的保持:\(\hat b_e=b_l/a_{le}\ge0\);对其余 \(i\),\(\hat b_i=b_i-a_{ie}(b_l/a_{le})\)——若 \(a_{ie}>0\),由最小比值 \(b_l/a_{le}\le b_i/a_{ie}\) 得 \(\hat b_i\ge0\);若 \(a_{ie}\le0\),\(\hat b_i\ge b_i\ge0\)。无界情形:所有 \(a_{ie}\le0\),令 \(x_e=t\to\infty\)、其余非基本变量为 0,基本变量 \(b_i-a_{ie}t\ge0\) 始终可行,目标值 \(\nu+c_et\to\infty\)。\(\square\)
最小比值检验是保持可行性的关键——选错出基变量会让某个 \(b_i\) 变负,基本解不再可行。
金融直觉:最小比值检验就是"找瓶颈"。你决定加大某个项目的投入 \(x_e\),每条约束都会说"我最多还能容忍你加到多少":第 \(i\) 条约束的剩余额度是 \(b_i\),每加一单位 \(x_e\) 消耗 \(a_{ie}\),所以它允许的上限是 \(b_i/a_{ie}\)。真正能加到的量是所有上限里最小的那个——最先被用完的资源决定了你能走多远,这条约束随即变"紧",它的松弛变量变成 0、出基。
如果 \(a_{ie}\le0\),说明加大 \(x_e\) 不消耗这项资源甚至还释放资源,它就不构成限制(比值记为 \(\infty\))。如果所有约束都不限制 \(x_e\),就可以无限加仓、目标无限增大——这就是"无界"。在组合优化里看到"unbounded",几乎总是漏写了满仓、杠杆或单票上限之类的约束。
29.4.5 退化与循环
每次转轴目标值不会下降(习题 29.3-2):\(\hat\nu=\nu+c_e\hat b_e\),\(c_e>0\),\(\hat b_e\ge0\)。目标值不变当且仅当 \(b_l=0\),这叫退化(degeneracy)。
退化例 \(z=x_1+x_2+x_3\),\(x_4=8-x_1-x_2\),\(x_5=x_2-x_3\)。\(x_1\) 入、\(x_4\) 出,得 \(z=8+x_3-x_4\),\(x_1=8-x_2-x_4\),\(x_5=x_2-x_3\)。此时只能 \(x_3\) 入、\(x_5\) 出,由于 \(b_5=0\),目标仍为 8:\(z=8+x_2-x_4-x_5\),\(x_3=x_2-x_5\)。再以 \(x_2\) 入、\(x_1\) 出,目标升到 16。
退化可能导致循环(cycling):若干次转轴后回到同一个松弛型,确定性的 SIMPLEX 就永远不会终止。
引理 29.4 给定基本变量集合 \(B\),松弛型唯一确定。(两个同 \(B\) 的松弛型相减,对每个 \(i\) 用代数引理 29.3——"若一个线性恒等式对所有取值成立,则对应系数相等、常数项为零"——即得。)
引理 29.5 \(n+m\) 个变量中选 \(m\) 个作基,至多 \(\binom{n+m}{m}\) 种。所以若 SIMPLEX 在 \(\binom{n+m}{m}\) 次迭代内不终止,它一定在循环。
防止循环:(a) 对 \(b\) 做微小扰动,使不同基本解目标值互不相同;(b) Bland 规则(Bland's rule):入基、出基出现平局时都选下标最小的变量。
引理 29.6、29.7 采用 Bland 规则时 SIMPLEX 必终止:要么报告无界,要么在至多 \(\binom{n+m}{m}\) 次迭代内返回可行解。
29.4.6 复杂度
每次转轴 \(O(mn)\);迭代次数最坏为 \(\binom{n+m}{m}\),可以是指数级。
29.4.7 最坏情况:Klee-Minty 立方体
Klee 与 Minty 构造了一族 LP,使得按"目标系数最大者入基"(Dantzig 规则)的单纯形法走遍立方体的全部 \(2^n\) 个顶点,需要 \(2^n-1\) 次转轴(原书注记)。29.8.1 节的代码复现了这一现象。实践中单纯形法很少表现出这种行为:Borgwardt 等证明了某些概率模型下的期望多项式时间;Spielman 与 Teng 的平滑分析(smoothed analysis)证明,对输入做微小随机扰动后,期望运行时间是多项式的,这解释了单纯形法的实际高效。
29.5 对偶性
29.5.1 动机
29.4 节证明了 SIMPLEX 会终止,却没有证明终止时的解最优。我们需要一个"证书"。第 26 章已经见过这样的证书:给定流 \(f\),若能找到容量等于 \(|f|\) 的割,\(f\) 就是最大流。线性规划对偶(linear-programming duality)把这个思想推广到一般 LP:给一个最大化问题(原始问题,primal),构造一个最小化问题(对偶问题,dual),使任何对偶可行解的目标值都是原始最优值的上界,并且两者最优值相等。
29.5.2 对偶的定义
标准型原始问题 \(\max c^Tx\),s.t. \(Ax\le b\),\(x\ge0\) 的对偶是
即 \(\min b^Ty\),s.t. \(A^Ty\ge c\),\(y\ge0\)。构造规则:max 变 min;右端常数与目标系数互换;\(\le\) 变 \(\ge\)。原始问题的每个约束对应一个对偶变量,对偶问题的每个约束对应一个原始变量。
怎样想到它:给第 \(i\) 个约束乘上非负权重 \(y_i\) 再相加,得 \(\sum_j(\sum_ia_{ij}y_i)x_j\le\sum_ib_iy_i\)。若每个 \(x_j\) 的系数都不小于 \(c_j\),那么(因为 \(x\ge0\))\(c^Tx\le\sum_j(\sum_ia_{ij}y_i)x_j\le b^Ty\)——\(b^Ty\) 就是原始目标的一个上界。对偶问题就是在寻找最紧的这种上界。
推导拆解:用下面的例题把这个想法算一遍。原始约束是 \(x_1+x_2+3x_3\le30\),\(2x_1+2x_2+5x_3\le24\),\(4x_1+x_2+2x_3\le36\),目标 \(3x_1+x_2+2x_3\)。
取权重 \(y=(0,\frac16,\frac23)\),即第一条不用,第二条乘 \(\frac16\),第三条乘 \(\frac23\),再相加。左边 \(x_1\) 的系数是 \(0+\frac26+\frac83=3\),\(x_2\) 的系数是 \(0+\frac26+\frac23=1\),\(x_3\) 的系数是 \(0+\frac56+\frac43=\frac{13}{6}\);右边是 \(0+4+24=28\)。于是 \(3x_1+x_2+\frac{13}{6}x_3\le28\)。
因为 \(x_3\ge0\) 且 \(2\le\frac{13}{6}\),所以 \(3x_1+x_2+2x_3\le3x_1+x_2+\frac{13}{6}x_3\le28\)。不用解任何方程,就证明了"目标不可能超过 28"。而 29.4.2 节已经找到目标恰为 28 的方案,所以 28 就是最优值。
这里要求 \(y_i\ge0\) 的原因:给"\(\le\)"不等式乘负数会把方向翻转,就不能再相加了。要求"每个 \(x_j\) 的系数 \(\ge c_j\)"的原因:只有这样,加权和才能"盖住"目标函数(这一步还用到了 \(x_j\ge0\))。对偶问题就是在所有合格的权重里挑出让右端最小、上界最紧的那一组。
例 29.4.2 节例题的对偶:
非标准型 LP 的对偶规则(习题 29.4-2,可由标准化推出)。对最大化的原始问题:
| 原始问题 | 对偶问题 |
|---|---|
| 第 \(i\) 个约束是 \(\le\) | \(y_i\ge0\) |
| 第 \(i\) 个约束是 \(=\) | \(y_i\) 无符号限制 |
| 第 \(i\) 个约束是 \(\ge\) | \(y_i\le0\) |
| \(x_j\ge0\) | 第 \(j\) 个对偶约束是 \(\ge c_j\) |
| \(x_j\) 无符号限制 | 第 \(j\) 个对偶约束是 \(=c_j\) |
| \(x_j\le0\) | 第 \(j\) 个对偶约束是 \(\le c_j\) |
对偶的对偶是原始问题(习题 29.4-5)。
29.5.3 弱对偶
引理 29.8(弱对偶,weak duality) 若 \(\bar x\) 原始可行、\(\bar y\) 对偶可行,则 \(\sum_jc_j\bar x_j\le\sum_ib_i\bar y_i\)。
证明:
第一步用对偶约束和 \(\bar x_j\ge0\),第二步用原始约束和 \(\bar y_i\ge0\)。\(\square\)
推论 29.9 若原始可行解与对偶可行解的目标值相等,则二者分别是原始、对偶的最优解。
第 26 章的推论 26.5"任意流值 ≤ 任意割容量"正是最大流问题的弱对偶(习题 29.4-6)。
29.5.4 从最终松弛型读出对偶解
例题的最终松弛型是 \(z=28-x_3/6-x_5/6-2x_6/3\)。规则:对偶变量等于最终目标函数中对应松弛变量系数的相反数:
例中 \(x_4\in B\),所以 \(\bar y_1=0\);\(\bar y_2=-c'_5=1/6\);\(\bar y_3=-c'_6=2/3\)。对偶目标值 \(30\cdot0+24\cdot\frac16+36\cdot\frac23=28\),与原始最优值相等。由推论 29.9,两者都是最优的——这就证明了 28 是最优值。
29.5.5 强对偶
定理 29.10(线性规划对偶定理) 设 SIMPLEX 对原始问题 \((A,b,c)\) 返回 \(\bar x\),\(N,B\) 为最终松弛型的非基本、基本变量集合,\(c'\) 为最终目标系数,\(\bar y\) 由 (29.91) 定义。则 \(\bar x\) 原始最优、\(\bar y\) 对偶最优,且 \(\sum_jc_j\bar x_j=\sum_ib_i\bar y_i\)。
证明:
- 最终目标 \(z=\nu'+\sum_{j\in N}c'_jx_j\),终止条件给出 \(c'_j\le0\)(\(j\in N\))。补定义 \(c'_j=0\)(\(j\in B\)),则 \(z=\nu'+\sum_{j=1}^{n+m}c'_jx_j\)。基本解中非基本变量为 0、基本变量系数为 0,所以原始目标值 \(\sum_jc_j\bar x_j=\nu'\)。
- 所有松弛型都等价,所以对任意 \(x\),\(\sum_{j=1}^nc_jx_j=\nu'+\sum_{j=1}^{n+m}c'_jx_j\)。把 \(x_{n+i}=b_i-\sum_ja_{ij}x_j\) 代入,并用 \(c'_{n+i}=-\bar y_i\),整理得
推导拆解:从第 2 步的恒等式到 (29.99) 只有一次代入和一次合并同类项。
把 \(\nu'+\sum_{j=1}^{n+m}c'_jx_j\) 拆成两部分:原变量部分 \(\sum_{j=1}^nc'_jx_j\),松弛变量部分 \(\sum_{i=1}^mc'_{n+i}x_{n+i}\)。
松弛变量部分用 \(c'_{n+i}=-\bar y_i\) 和 \(x_{n+i}=b_i-\sum_ja_{ij}x_j\) 代入: \(\sum_ic'_{n+i}x_{n+i}=\sum_i(-\bar y_i)\big(b_i-\sum_ja_{ij}x_j\big)=-\sum_ib_i\bar y_i+\sum_j\big(\sum_ia_{ij}\bar y_i\big)x_j\)。 最后一步交换了两个求和的顺序(有限项求和,先按 \(i\) 加还是先按 \(j\) 加结果相同)。
把常数项归到一起(\(\nu'-\sum_ib_i\bar y_i\)),把 \(x_j\) 的系数归到一起(\(c'_j+\sum_ia_{ij}\bar y_i\)),就是 (29.99)。
第 3 步"恒等式两边系数相等"的意思是:这个等式对任意 \(x\) 都成立,那么令 \(x=0\) 可知常数项为 0;再令 \(x\) 只有第 \(j\) 个分量为 1,可知 \(x_j\) 的系数相等。这和"两个多项式处处相等,则各项系数相等"是同一个道理。
- 这是关于 \(x_1..x_n\) 的恒等式,由引理 29.3,常数项 \(\nu'-\sum_ib_i\bar y_i=0\),且 \(c'_j+\sum_ia_{ij}\bar y_i=c_j\)。前者说明对偶目标值等于原始最优值 \(\nu'\)。
- 对偶可行性:\(c'_j\le0\) 对所有 \(j\) 成立,所以 \(c_j=c'_j+\sum_ia_{ij}\bar y_i\le\sum_ia_{ij}\bar y_i\);又 \(\bar y_i=-c'_{n+i}\ge0\)。由推论 29.9 得最优。\(\square\)
这个证明是构造性的:单纯形法在求出原始最优解的同时,免费给出了对偶最优解。
29.5.6 影子价格
对偶变量有清楚的经济含义:\(\bar y_i\) 是第 \(i\) 种"资源"(约束右端 \(b_i\))每增加一单位,最优目标值的增加量,称为影子价格(shadow price)。例题中把 \(b_2\) 从 24 增到 25,最优值从 28 变为 \(28+1/6=169/6\);把 \(b_3\) 从 36 增到 37,最优值变为 \(28+2/3=86/3\)(29.8 节已用代码核对);而 \(b_1\) 有 18 的富余,再增加也没有用,\(\bar y_1=0\)。
两点限定:影子价格只在最优基不变的范围内有效(右端变化太大会换基,斜率随之改变,最优值是右端的分段线性凹函数);在退化最优解处,对偶解可能不唯一,"左导数"和"右导数"不同。
量化解读:组合优化中,对偶变量告诉你每条约束(行业中性、单票上限、换手上限、跟踪误差预算)"放松一单位能多挣多少目标值"。这是约束归因(constraint attribution)的基础:一条影子价格很高的约束在大量消耗 alpha,值得重新审视它是否必要。
金融直觉:回到资本限额下的资本预算。设变量是各项目的参与比例,约束是"第 1 年资本 ≤ 1 亿""第 2 年资本 ≤ 8000 万""工程团队工时 ≤ 5 万小时",目标是总 NPV。求解后,每条约束都有一个对偶变量:
- 第 1 年资本的影子价格若是 0.18,表示多 1 元预算能多 0.18 元 NPV。只要每追加 1 元融资带来的额外代价(按现值计,例如发行费用、超出原资本成本的那部分利息)低于 0.18 元,融资就值得。
- 工时约束的影子价格若是 0,说明工时有富余(互补松弛:有余量的资源边际价值为零),招人不能增加价值。
- 影子价格还给出项目的"内部定价":一个未入选项目的检验数 \(c_j-\sum_ia_{ij}\bar y_i\),就是"它的 NPV 减去它按影子价格计算的资源占用成本"。为负就不该投,这是盈利指数法在多约束下的正确推广。
原文的"两点限定"也很重要:影子价格是局部的边际量。预算从 1 亿加到 1.01 亿可以用它估计,从 1 亿加到 2 亿就不行,因为中途最优项目组合会改变,边际价值也随之下降(最优值是右端的凹函数,和"边际收益递减"一致)。
29.5.7 互补松弛
思考题 29-2(互补松弛,complementary slackness) 设 \(\bar x,\bar y\) 分别原始可行、对偶可行,则二者同为最优当且仅当
直观地说:有富余的资源,影子价格为零;影子价格为正的资源,一定被用光。证明只需注意弱对偶证明中的两个不等号同时取等。例题中第一个约束有 18 的松弛,\(\bar y_1=0\);\(x_3=0\),第三个对偶约束 \(3\bar y_1+5\bar y_2+2\bar y_3=13/6>2\) 有 \(1/6\) 的松弛。
互补松弛是检验求解器输出、以及从原始解反推对偶解的实用工具;它也是内点法和 KKT 条件(第 04 册第 12 章)在线性情形下的形式。
29.6 初始基本可行解
29.6.1 问题
LP 可行,但令所有原变量为 0 的初始基本解不一定可行。例:
\(x_1=x_2=0\) 违反第二个约束(\(b_2<0\))。
29.6.2 辅助线性规划
引理 29.11 设 \(L\) 为标准型 LP,引入新变量 \(x_0\),定义
则 \(L\) 可行当且仅当 \(L_{aux}\) 的最优值为 0。
证明:\(L\) 的可行解配上 \(x_0=0\) 是 \(L_{aux}\) 的可行解,目标值 0,而 \(-x_0\le0\),所以最优。反之,最优值为 0 意味着 \(\bar x_0=0\),其余分量满足 \(L\) 的约束。\(\square\)
直观上,\(x_0\) 是"允许每个约束被违反的量",最小化它就是在问"能不能一点都不违反"。
白话解释:为什么需要这一节?单纯形法要求从一个"角"出发,最省事的角是"所有决策变量取 0"。但如果某条约束是"至少要拿到 2.5 万张农村选票"或"至少持有 30% 债券",全取 0 就违反了约束,起点本身不可行。辅助问题的做法是:先给每条约束同样宽限 \(x_0\),宽限足够大时全取 0 一定可行;然后用单纯形法去最小化宽限。宽限能压到 0,说明原问题有可行点,而且此时恰好停在原问题的一个角上,可以接着做第二阶段;压不到 0,说明约束之间互相矛盾,原问题不可行。实务中求解器报告 "infeasible" 时,走的就是这个逻辑。
29.6.3 INITIALIZE-SIMPLEX
def INITIALIZE_SIMPLEX(A, b, c):
k = argmin_i b[i]
if b[k] >= 0: # 初始基本解已可行
return ({1..n}, {n+1..n+m}, A, b, c, 0)
构造 L_aux 的松弛型:非基 {0,1..n},基 {n+1..n+m}
l = n + k # 基本解中最负的那个基变量出基
PIVOT(..., l, 0) # x0 入基,之后基本解可行
反复执行 SIMPLEX 的 while 循环,直到求得 L_aux 的最优解
if 最优解中 x0 == 0:
if x0 是基变量:
做一次退化转轴,选任一 a[0][e] != 0 的 e 入基,使 x0 出基
从约束中删去 x0;恢复 L 的原目标函数,
并把其中出现的每个基变量替换为其约束的右端表达式
return 修改后的最终松弛型
else:
return "infeasible"
关键技巧:只需一次转轴(\(x_0\) 入基、最负的 \(x_{n+k}\) 出基)就能让 \(L_{aux}\) 的基本解可行。因为 \(L_{aux}\) 中 \(x_0\) 的系数 \(a_{i0}=-1\),转轴后 \(\bar x_0=-b_k>0\),其余 \(\bar x_i=b_i-b_k\ge0\)(\(b_k\) 是最小值)。
例 辅助问题的松弛型 \(z=-x_0\),\(x_3=2-2x_1+x_2+x_0\),\(x_4=-4-x_1+5x_2+x_0\),基本解 \(x_4=-4\) 不可行。
- \(x_0\) 入基、\(x_4\) 出基:\(z=-4-x_1+5x_2-x_4\),\(x_0=4+x_1-5x_2+x_4\),\(x_3=6-x_1-4x_2+x_4\);基本解 \((x_0,\dots,x_4)=(4,0,0,6,0)\) 可行。
- \(x_2\) 入基、\(x_0\) 出基:\(z=-x_0\),\(x_2=\frac45-\frac{x_0}{5}+\frac{x_1}{5}+\frac{x_4}{5}\),\(x_3=\frac{14}{5}+\frac{4x_0}{5}-\frac{9x_1}{5}+\frac{x_4}{5}\)。最优值 0,原问题可行。
- 删去 \(x_0\),恢复原目标 \(2x_1-x_2\) 并代入 \(x_2\) 的表达式:\(z=-\frac45+\frac{9x_1}{5}-\frac{x_4}{5}\),约束 \(x_2=\frac45+\frac{x_1}{5}+\frac{x_4}{5}\),\(x_3=\frac{14}{5}-\frac{9x_1}{5}+\frac{x_4}{5}\)。这个松弛型的基本解可行,交给 SIMPLEX 继续求解,最终得 \(x=(14/9,10/9)\),目标值 2(29.8.1 节代码输出)。
引理 29.12 若 \(L\) 不可行,INITIALIZE-SIMPLEX 返回"infeasible";否则返回一个基本解可行的合法松弛型。(不可行时 \(L_{aux}\) 的最优值为负且有限;在 INITIALIZE-SIMPLEX 内运行主循环也不会返回"unbounded",因为 \(-x_0\le0\) 有上界,习题 29.5-2。)
这就是两阶段法(two-phase method):第一阶段求可行基,第二阶段求最优。
29.6.4 线性规划基本定理
定理 29.13(Fundamental theorem of linear programming) 标准型 LP \(L\) 必为以下三者之一:
- 有有限最优值的最优解;
- 不可行;
- 无界。
SIMPLEX 在三种情形下分别返回最优解、"infeasible"、"unbounded"。
一个推论值得记住:不存在"最优值有限但取不到"的 LP——这与非线性优化不同,也是严格不等式被排除的原因。
29.7 思考题与注记
- 29-1 线性不等式可行性:(a) 用 LP 算法解可行性问题(目标取 0);(b) 反过来用可行性算法解 LP——把原始约束、对偶约束和"原始目标 ≥ 对偶目标"放在一起,由强对偶,它的任何可行解都给出最优解。两者在多项式意义下等价。
- 29-2 互补松弛:见 29.5.7 节。
- 29-3 整数线性规划:判定可行性就是 NP 难的。(a) 弱对偶仍成立;(b) 强对偶不一定成立;(c) 记整数原始、LP 松弛原始、LP 松弛对偶、整数对偶的最优值为 \(IP,P,D,ID\),则 \(IP\le P=D\le ID\)。\(P-IP\) 叫整数间隙(integrality gap)。这解释了为何"先解 LP 松弛再取整"可能明显次优。
- 29-4 Farkas 引理:对 \(m\times n\) 矩阵 \(A\) 和 \(n\) 维向量 \(c\),以下两个系统恰有一个有解:
这是对偶理论的"择一定理"形式。它的金融版本就是资产定价基本定理(29.8.3 节)。
- 29-5 最小费用循环流:没有源汇和需求,只有容量约束和处处守恒,求费用最小的流。所有边费用为正时最优解是零流;最大流可化为最小费用循环流(加一条 \(t\to s\) 的大容量负费用边);单源最短路也可以。它是网络优化问题的统一框架。
注记:单纯形法由 G. Dantzig 于 1947 年提出,至今其变体仍是最常用的 LP 方法。第一个多项式时间算法是 Khachian(1979)的椭球法;Karmarkar 提出第一个内点法。Klee 与 Minty 构造了需要 \(2^n-1\) 次迭代的例子;Spielman 与 Teng 用平滑分析解释了其实际高效。网络单纯形法(network simplex)在最短路、最大流、最小费用流上有多项式时间的变体。原书的讲法取自 Chvátal 的教材。
29.8 量化实战
29.8.1 单纯形法的实现与验证
下面是原书算法的直译实现:PIVOT、带 Bland 规则(或 Dantzig 最大系数规则)的主循环、两阶段 INITIALIZE-SIMPLEX,以及按式 (29.91) 读出对偶解。用 Fraction 做精确有理数运算,这样能和原书的分数逐一对照。
"""原书第 29 章单纯形算法的直译实现(松弛型 + PIVOT + 两阶段初始化 + 对偶解)。
变量编号沿用原书:原变量 1..n,松弛变量 n+1..n+m,辅助变量 0。用 Fraction 做精确运算。"""
from fractions import Fraction as Fr
def pivot(N, B, A, b, c, v, l, e):
"""原书 PIVOT:x_e 入基,x_l 出基。A[i][j] 是松弛型 x_i = b_i - sum A[i][j] x_j 中的系数。"""
An, bn, cn = {}, {}, {}
bn[e] = b[l] / A[l][e]
An[e] = {j: A[l][j] / A[l][e] for j in N if j != e}
An[e][l] = 1 / A[l][e]
for i in B:
if i == l: continue
bn[i] = b[i] - A[i][e] * bn[e]
An[i] = {j: A[i][j] - A[i][e] * An[e][j] for j in N if j != e}
An[i][l] = -A[i][e] * An[e][l]
vn = v + c[e] * bn[e]
cn = {j: c[j] - c[e] * An[e][j] for j in N if j != e}
cn[l] = -c[e] * An[e][l]
return (N - {e}) | {l}, (B - {l}) | {e}, An, bn, cn, vn
def _loop(N, B, A, b, c, v, rule="bland", trace=None):
"""SIMPLEX 的 while 循环。rule='bland':最小下标;'dantzig':目标系数最大者入基。"""
pivots = 0
while True:
cand = sorted(j for j in N if c[j] > 0)
if not cand:
return N, B, A, b, c, v, pivots
e = cand[0] if rule == "bland" else max(cand, key=lambda j: (c[j], -j))
ratios = [(b[i] / A[i][e], i) for i in sorted(B) if A[i][e] > 0] # 最小比值检验
if not ratios:
raise ArithmeticError("unbounded")
l = min(ratios)[1] # 平局取最小下标
N, B, A, b, c, v = pivot(N, B, A, b, c, v, l, e)
pivots += 1
if trace is not None:
trace.append((e, l, v))
def initialize_simplex(A0, b0, c0):
m, n = len(b0), len(c0)
N, B = set(range(1, n + 1)), set(range(n + 1, n + m + 1))
A = {n + 1 + i: {j + 1: Fr(A0[i][j]) for j in range(n)} for i in range(m)}
b = {n + 1 + i: Fr(b0[i]) for i in range(m)}
c = {j + 1: Fr(c0[j]) for j in range(n)}
k = min(b, key=lambda i: (b[i], i))
if b[k] >= 0:
return N, B, A, b, c, Fr(0)
# 辅助线性规划 L_aux:每个约束加 -x0,目标 max -x0
Na = N | {0}
for i in B: A[i][0] = Fr(-1)
ca = {j: Fr(0) for j in N}; ca[0] = Fr(-1)
Na, Ba, A, b, ca, va = pivot(Na, B, A, b, ca, Fr(0), k, 0) # 一次转轴即得可行基本解
Na, Ba, A, b, ca, va, _ = _loop(Na, Ba, A, b, ca, va)
if va != 0: # 最优值 < 0:原问题不可行
raise ArithmeticError("infeasible")
if 0 in Ba: # 退化:x0 仍在基中,做一次退化转轴
e = min(j for j in Na if A[0][j] != 0)
Na, Ba, A, b, ca, va = pivot(Na, Ba, A, b, ca, va, 0, e)
Na = Na - {0}
for i in Ba: A[i].pop(0, None)
# 恢复原目标函数,并把其中的基变量用其约束右端代入
cn, v = {j: Fr(0) for j in Na}, Fr(0)
for j in range(1, n + 1):
cj = Fr(c0[j - 1])
if j in Na:
cn[j] += cj
else:
v += cj * b[j]
for k2 in Na: cn[k2] -= cj * A[j][k2]
return Na, Ba, A, b, cn, v
def simplex(A0, b0, c0, rule="bland", trace=None):
"""max c^T x s.t. A x <= b, x >= 0。返回 (x, 目标值, 对偶解 y, 转轴次数)。"""
m, n = len(b0), len(c0)
N, B, A, b, c, v = initialize_simplex(A0, b0, c0)
N, B, A, b, c, v, piv = _loop(N, B, A, b, c, v, rule, trace)
x = [b[j] if j in B else Fr(0) for j in range(1, n + 1)]
y = [-c[n + i] if (n + i) in N else Fr(0) for i in range(1, m + 1)] # 式 (29.91)
return x, v, y, piv
把上面的代码存为 simplex.py,再运行下面的验证脚本:
import numpy as np
from scipy.optimize import linprog
from simplex import simplex
fmt = lambda xs: "(" + ", ".join(str(x) for x in xs) + ")"
# 29.3 节例题 (29.53)-(29.57):max 3x1 + x2 + 2x3
A = [[1, 1, 3], [2, 2, 5], [4, 1, 2]]; b = [30, 24, 36]; c = [3, 1, 2]
for rule in ("dantzig", "bland"): # 原书例题的转轴顺序 = 最大系数规则
trace = []
x, z, y, _ = simplex(A, b, c, rule=rule, trace=trace)
print(rule, ":", ";".join(f"x{e} 入基 x{l} 出基 -> {v}" for e, l, v in trace))
print("最优解 x =", fmt(x), " z =", z, " 对偶解 y =", fmt(y), " 对偶目标 =", sum(bi * yi for bi, yi in zip(b, y)))
# 互补松弛(思考题 29-2):原约束有松弛 <=> 对应 y_i = 0;x_j > 0 <=> 对偶约束取等
slack = [bi - sum(a * xi for a, xi in zip(row, x)) for row, bi in zip(A, b)]
dual_slack = [sum(A[i][j] * y[i] for i in range(3)) - c[j] for j in range(3)]
print("原约束松弛量", fmt(slack), " 对偶约束松弛量", fmt(dual_slack))
# 29.5 节例题:初始基本解不可行,需要两阶段
x, z, y, _ = simplex([[2, -1], [1, -5]], [2, -4], [2, -1])
print("29.5 节例题:x =", fmt(x), " z =", z)
# 不可行与无界
for name, (A_, b_, c_) in {"习题29.1-6": ([[1, 1], [-2, -2]], [2, -10], [1, -2]),
"无界": ([[1, -1]], [1], [1, 1])}.items():
try:
simplex(A_, b_, c_)
except ArithmeticError as err:
print(name, "->", err)
# 政治问题 (29.6)-(29.10):min 广告费,>= 约束取负化为标准型
P = [[-2, 8, 0, 10], [5, 2, 0, 0], [3, -5, 10, -2]]; need = [50, 100, 25]
A_std = [[-a for a in row] for row in P]; b_std = [-v for v in need]
x, z, y, _ = simplex(A_std, b_std, [-1, -1, -1, -1])
print("政治问题:x =", fmt(x), f" 最少花费 = {-z} = {float(-z):.4f} 千美元,影子价格 y =", fmt(y))
res = linprog([1, 1, 1, 1], A_ub=A_std, b_ub=b_std, method="highs")
print(f" HiGHS 核对:{res.fun:.4f},约束对偶 {(-res.ineqlin.marginals).round(4)}")
# Klee-Minty 立方体:按"最大系数"规则选入基变量需要 2^n - 1 次转轴
def klee_minty(n):
A = [[2 ** (i - j + 1) if j < i else (1 if j == i else 0) for j in range(n)] for i in range(n)]
return A, [5 ** (i + 1) for i in range(n)], [2 ** (n - 1 - j) for j in range(n)]
for n in (3, 5, 8, 10):
A_, b_, c_ = klee_minty(n)
_, zd, _, pd = simplex(A_, b_, c_, rule="dantzig")
_, zb, _, pb = simplex(A_, b_, c_, rule="bland")
print(f"Klee-Minty n={n:2d}: 最大系数规则 {pd:4d} 次转轴(2^n-1={2**n-1}),Bland 规则 {pb:4d} 次,最优值 {zd}")
输出:
dantzig : x1 入基 x6 出基 -> 27;x3 入基 x5 出基 -> 111/4;x2 入基 x3 出基 -> 28
bland : x1 入基 x6 出基 -> 27;x2 入基 x5 出基 -> 28
最优解 x = (8, 4, 0) z = 28 对偶解 y = (0, 1/6, 2/3) 对偶目标 = 28
原约束松弛量 (18, 0, 0) 对偶约束松弛量 (0, 0, 1/6)
29.5 节例题:x = (14/9, 10/9) z = 2
习题29.1-6 -> infeasible
无界 -> unbounded
政治问题:x = (2050/111, 425/111, 0, 625/111) 最少花费 = 3100/111 = 27.9279 千美元,影子价格 y = (25/222, 23/111, 7/111)
HiGHS 核对:27.9279,约束对偶 [0.1126 0.2072 0.0631]
Klee-Minty n= 3: 最大系数规则 7 次转轴(2^n-1=7),Bland 规则 5 次,最优值 125
Klee-Minty n= 5: 最大系数规则 31 次转轴(2^n-1=31),Bland 规则 15 次,最优值 3125
Klee-Minty n= 8: 最大系数规则 255 次转轴(2^n-1=255),Bland 规则 67 次,最优值 390625
Klee-Minty n=10: 最大系数规则 1023 次转轴(2^n-1=1023),Bland 规则 177 次,最优值 9765625
读结果:
- 用最大系数规则,转轴序列 \(27\to111/4\to28\) 与原书逐步一致;用 Bland 规则,第二步选的是下标更小的 \(x_2\),两步就到最优。不同的入基规则走的是可行域上不同的顶点路径,终点相同。
- 对偶解 \((0,1/6,2/3)\)、对偶目标 28 与 29.5.4 节一致;互补松弛条件成立:第一个原始约束有 18 的松弛,对应 \(y_1=0\);\(x_3=0\),对应的第三个对偶约束有 \(1/6\) 的松弛。
- 29.6 节的两阶段例题得到 \(x=(14/9,10/9)\)、目标值 2;不可行(习题 29.1-6 的约束 \(x_1+x_2\le2\) 与 \(x_1+x_2\ge5\) 矛盾,目标任取)和无界的情形都被正确识别。
- 政治问题的最少花费是 \(3100/111\approx27.93\) 千美元,三个选区约束都是紧的。影子价格 \((25/222,23/111,7/111)\) 表示:郊区的选票要求每提高 1 千张,最少花费增加约 0.207 千美元,是三者中最"贵"的。HiGHS 求解器给出的对偶完全一致。
- Klee-Minty 立方体上,最大系数规则恰好需要 \(2^n-1\) 次转轴;Bland 规则在这个例子上少得多,但它同样有指数级的最坏例子。工业求解器使用更精细的定价规则(如最陡边规则)和扰动技术。
29.8.2 部分复制的指数跟踪:L1 跟踪误差的 LP
问题。某指数有 300 只成分股,但出于交易成本和管理成本,只想持有市值最大的 80 只,构造一个尽量贴近指数的组合。这叫部分复制(partial replication)或抽样复制。常见约束:满仓、不卖空、单票上限、行业权重偏离不超过 ±0.5%。
为什么用 L1。跟踪误差通常用主动收益的标准差(L2)度量,最小化它是一个二次规划(第 04 册第 16a、16b 章)。改用平均绝对偏差(mean absolute deviation,MAD):
就能写成 LP(这一思路即 Konno 与 Yamazaki 的 MAD 组合模型)。技巧是把每个绝对值拆成两个非负变量:
在最优解处,\(u_t\) 和 \(v_t\) 至多一个为正(否则同减一个正数,目标更小),所以 \(u_t+v_t=|\cdot|\)。L1 对极端日的权重比 L2 小,结果更稳健;代价是没有 L2 那样的解析结构。
import numpy as np
from scipy.optimize import linprog
from scipy.sparse import hstack, vstack, csr_matrix, identity
rng = np.random.default_rng(29)
# ---------- 模拟市场:300 只股票、10 个行业,单因子 + 行业因子 + 特质收益 ----------
N, S, T = 300, 10, 1250
sector = np.repeat(np.arange(S), N // S)
beta = rng.uniform(0.7, 1.3, N)
f_m = rng.normal(0.0003, 0.010, T); f_s = rng.normal(0, 0.006, (T, S))
idio = rng.normal(0, 1, (T, N)) * rng.uniform(0.010, 0.022, N)
R = f_m[:, None] * beta + f_s[:, sector] + idio
cap = rng.lognormal(0, 1.2, N); w_idx = cap / cap.sum()
r_idx = R @ w_idx
tr, te = slice(0, 1000), slice(1000, 1250) # 前 1000 天估计,后 250 天检验
elig = np.argsort(-cap)[:80] # 只允许持有市值最大的 80 只(部分复制)
n = len(elig); sec_e = sector[elig]
W_sec = np.bincount(sector, weights=w_idx, minlength=S) # 指数的行业权重
def tracking_lp(alpha=None, mad_budget=None, sec_tol=0.005, ub=0.04):
"""变量 [w (n), u (T), v (T)]。R w - r_idx = u - v,MAD = mean(u + v)。
alpha 为 None 时最小化 MAD;否则在 MAD <= mad_budget 下最大化 alpha'w。"""
Rt, rt = R[tr][:, elig], r_idx[tr]; Tt = len(rt)
I = identity(Tt, format="csr")
A_eq = vstack([hstack([csr_matrix(Rt), -I, I]),
csr_matrix(np.r_[np.ones(n), np.zeros(2 * Tt)])])
b_eq = np.r_[rt, 1.0]
G = np.zeros((S, n)); G[sec_e, np.arange(n)] = 1 # 行业权重 = G w
Z = np.zeros((S, 2 * Tt))
A_ub = [np.hstack([G, Z]), np.hstack([-G, Z])] # |G w - W_sec| <= sec_tol
b_ub = [W_sec + sec_tol, -(W_sec - sec_tol)]
mad_row = np.r_[np.zeros(n), np.ones(2 * Tt) / Tt]
if alpha is None:
cost = mad_row
else:
cost = np.r_[-alpha, np.zeros(2 * Tt)]
A_ub.append(mad_row[None, :]); b_ub.append([mad_budget])
bounds = [(0, ub)] * n + [(0, None)] * (2 * Tt)
res = linprog(cost, A_ub=np.vstack(A_ub), b_ub=np.concatenate(b_ub),
A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="highs")
assert res.status == 0, res.message
return res
def te_stats(w, sl):
active = R[sl][:, elig] @ w - r_idx[sl]
return np.abs(active).mean() * 1e4, active.std() * np.sqrt(252) * 100
res = tracking_lp()
w = res.x[:n]
naive = w_idx[elig] / w_idx[elig].sum() # 对照:前 80 只按市值重新归一
print(f"80 只成分股覆盖指数权重 {w_idx[elig].sum():.1%};LP 规模:{res.x.size} 个变量,HiGHS 迭代 {res.nit} 次")
for name, ww in (("市值归一", naive), ("LP 最小 MAD", w)):
(m1, t1), (m2, t2) = te_stats(ww, tr), te_stats(ww, te)
print(f"{name:8s}: 样本内 MAD {m1:5.2f}bp / 跟踪误差 {t1:4.2f}%;样本外 MAD {m2:5.2f}bp / 跟踪误差 {t2:4.2f}%")
print(f"LP 持仓数 {np.sum(w > 1e-6)},触及 4% 上限 {np.sum(w > 0.04 - 1e-9)} 只")
# 影子价格:marginals = d(最优目标)/d(约束右端)
m_up, m_lo = res.ineqlin.marginals[:S], res.ineqlin.marginals[S:2 * S]
for s in np.nonzero((np.abs(m_up) > 1e-10) | (np.abs(m_lo) > 1e-10))[0]:
side = "上限" if abs(m_up[s]) > 1e-10 else "下限"
mu = m_up[s] if side == "上限" else m_lo[s] # 两种约束都写成 <=,放宽即右端增大
print(f" 行业 {s} 的{side}约束是紧的:放宽 1 个百分点,MAD 变化 {mu * 0.01 * 1e4:+.3f}bp")
ubm = res.upper.marginals[:n]
for i in np.nonzero(ubm < -1e-10)[0]:
print(f" 股票 {elig[i]} 的 4% 上限是紧的:放宽 1 个百分点,MAD 变化 {ubm[i] * 0.01 * 1e4:+.3f}bp")
# ---------- 增强指数:MAD 预算下最大化 alpha,对偶 = 前沿斜率 ----------
alpha = rng.normal(0, 0.02, n) # 年化 alpha 预测
mad0 = res.fun
print("\nMAD 预算(bp) 组合alpha(%) 预算约束影子价格 有限差分斜率")
for k in (1.05, 1.2, 1.5, 2.0):
tau = mad0 * k
r1 = tracking_lp(alpha, tau); r2 = tracking_lp(alpha, tau * 1.001)
shadow = -r1.ineqlin.marginals[-1] # linprog 最小化 -alpha'w,取负还原
fd = (-r2.fun + r1.fun) / (tau * 0.001)
print(f"{tau * 1e4:10.3f} {-r1.fun * 100:11.3f} {shadow:16.2f} {fd:14.2f}")
# ---------- 状态价格与无套利(Farkas 引理 / 思考题 29-4) ----------
D = np.array([[1, 1, 1], # 债券
[80, 100, 120], # 股票
[0, 10, 30], # 行权价 90 的看涨期权
[0, 0, 10]]) # 行权价 110 的看涨期权
for p in (np.array([0.95, 94, 11.5, 2.5]), np.array([0.95, 94, 10.0, 2.5])): # 第二组:C90 报价偏低
st = linprog(np.zeros(3), A_eq=D, b_eq=p, bounds=[(0, None)] * 3, method="highs")
if st.status == 0:
print("价格", p, "-> 存在状态价格 psi =", st.x.round(4), ",无套利")
else: # Farkas:找组合 theta,使到期在每个状态都不亏(D^T theta >= 0),而今天成本 p'theta < 0
# 固定买入 1 份 C90,其余头寸自由,求今天成本最低且到期各状态都不亏的组合
arb = linprog(p, A_ub=-D.T, b_ub=np.zeros(3), bounds=[(None, None)] * 2 + [(1, 1), (None, None)],
method="highs")
print("价格", p, "-> 无状态价格;套利组合 theta =", arb.x.round(4),
f",今天净收入 {-arb.fun:.3f},各状态到期收益 {(D.T @ arb.x).round(6) + 0}")
输出:
80 只成分股覆盖指数权重 73.9%;LP 规模:2080 个变量,HiGHS 迭代 1846 次
市值归一 : 样本内 MAD 6.36bp / 跟踪误差 1.26%;样本外 MAD 6.27bp / 跟踪误差 1.23%
LP 最小 MAD: 样本内 MAD 4.94bp / 跟踪误差 1.04%;样本外 MAD 5.76bp / 跟踪误差 1.12%
LP 持仓数 80,触及 4% 上限 2 只
行业 1 的下限约束是紧的:放宽 1 个百分点,MAD 变化 -0.017bp
行业 2 的上限约束是紧的:放宽 1 个百分点,MAD 变化 -0.036bp
行业 3 的上限约束是紧的:放宽 1 个百分点,MAD 变化 -0.022bp
行业 7 的下限约束是紧的:放宽 1 个百分点,MAD 变化 -0.001bp
行业 9 的下限约束是紧的:放宽 1 个百分点,MAD 变化 -0.045bp
股票 258 的 4% 上限是紧的:放宽 1 个百分点,MAD 变化 -0.695bp
股票 286 的 4% 上限是紧的:放宽 1 个百分点,MAD 变化 -0.253bp
MAD 预算(bp) 组合alpha(%) 预算约束影子价格 有限差分斜率
5.190 0.329 50.05 50.05
5.931 0.581 27.15 27.15
7.414 0.904 18.26 18.26
9.885 1.275 12.23 12.23
价格 [ 0.95 94. 11.5 2.5 ] -> 存在状态价格 psi = [0.3 0.4 0.25] ,无套利
价格 [ 0.95 94. 10. 2.5 ] -> 无状态价格;套利组合 theta = [40. -0.5 1. -1. ] ,今天净收入 1.500,各状态到期收益 [0. 0. 0.]
读结果:跟踪组合
- 样本内外。按市值把 80 只股票重新归一,样本外日均绝对偏差 6.27 bp、年化跟踪误差 1.23%;LP 组合在样本内把 MAD 降到 4.94 bp,样本外 5.76 bp、1.12%。样本外的改善小于样本内——LP 在 1000 天里用 80 个权重拟合了一部分噪声,这是所有"用历史收益优化权重"的方法都有的估计误差问题。训练窗口越短、可调权重越多,过拟合越严重(读者可以把训练期改成 250 天,会看到样本外反而不如市值归一)。实务中常加的正则化包括:让权重贴近市值权重、用因子模型的协方差代替原始收益样本、限制换手。
- 影子价格。HiGHS 的
marginals就是最优目标对约束右端的导数(式 29.91 的数值版本)。五个行业约束是紧的:例如行业 9 的下限每放宽 1 个百分点,MAD 降低 0.045 bp;两只股票的 4% 上限也是紧的,其中股票 258 的上限最"贵",放宽 1 个百分点可降低 MAD 0.695 bp——它在指数中的权重高于 4%,上限迫使组合低配它。这类数字就是约束归因:哪些约束在消耗跟踪精度,值不值得为它付出代价。
读结果:增强指数前沿
把问题倒过来:在 MAD 不超过预算 \(\tau\) 的前提下最大化预期 alpha。MAD 预算这条约束的影子价格,就是 alpha-跟踪误差前沿的斜率 \(d\alpha/d\tau\)。表中影子价格与有限差分斜率吻合到小数点后两位,验证了 29.5.6 节的解释。斜率从 50 降到 12:最初每多 1 bp 的 MAD 预算能换来约 0.5% 的年化 alpha,之后越来越少——LP 的最优值是右端的凹分段线性函数,这就是"边际效用递减"在 LP 里的精确形式。基金经理据此可以判断:在哪个跟踪误差水平上,再放宽预算已经不划算。
29.8.3 Farkas 引理与无套利
单期市场有 \(S\) 个状态、\(n\) 个资产,资产 \(i\) 在状态 \(s\) 的到期收益为 \(D_{is}\),今天价格为 \(p_i\)。状态价格(state price)向量 \(\psi\ge0\) 若满足 \(D\psi=p\),就给出了一个与所有报价一致的线性定价规则。Farkas 引理(思考题 29-4)说,以下两者恰有一个成立:
(I) 就是套利:一个今天净收入为正(成本为负)、到期在每个状态都不亏的组合 \(\theta\)。于是"不存在这种套利 ⟺ 存在非负的状态价格"——资产定价基本定理的离散版本(第 08 册)。(严格的版本要求 \(\psi>0\) 以排除"不花钱、可能赚钱"的弱套利,证明思路相同。)
金融直觉:这和 CFA 里的单步二叉树是同一件事。二叉树有两个状态(上涨、下跌),状态价格 \(\psi_u,\psi_d\) 满足"债券价格 \(=\psi_u+\psi_d\)""股票价格 \(=\psi_uS_u+\psi_dS_d\)"。两个方程、两个未知数,正好解出唯一的 \(\psi\);再除以 \(\psi_u+\psi_d\) 就得到风险中性概率 \(q=\psi_u/(\psi_u+\psi_d)\)。CFA 要求 \(S_d<S_0(1+r)<S_u\),正是为了保证解出来的 \(\psi_u,\psi_d\) 都为正,即无套利。
Farkas 引理把这个结论推广到任意多个状态、任意多个资产:方程组 \(D\psi=p\) 未必有唯一解(不完全市场)、也可能根本没有非负解。它保证两种情况必居其一,没有第三种:要么能找到一组非负的状态价格把所有报价同时"解释"出来,要么能找到一个套利组合。LP 求解器在前一种情况下给出 \(\psi\),在后一种情况下(通过对偶)给出套利组合 \(\theta\),所以代码只需调用一次
linprog就能判断一组报价是否自洽。
代码中第一组价格存在状态价格 \(\psi=(0.3,0.4,0.25)\),无套利。第二组把行权价 90 的看涨期权报价从 11.5 压到 10:债券、股票、行权价 110 的期权三者已经唯一确定了 \(\psi\)(完全市场),C90 的价格必须是 \(0\cdot0.3+10\cdot0.4+30\cdot0.25=11.5\)。LP 找到的套利组合是:买入 1 份 C90、卖出 1 份 C110、卖空 0.5 股股票、买入 40 份债券。到期在三个状态的收益都是 0,今天却净收入 1.5——正好是错价的幅度。这个组合就是"用其他资产复制 C90 并反向做"的静态套利,LP 自动把它找了出来。同样的方法可以检验一张期权链的报价是否存在静态套利,或者在不完全市场中用 LP 求某个衍生品的无套利价格上下界(超复制问题及其对偶)。
29.8.4 实务要点
- 用求解器,但要懂它的输出。HiGHS(
scipy.optimize.linprog的默认后端)、Gurobi、CPLEX 都同时提供单纯形法和内点法。读懂"不可行""无界""对偶值""基状态"这些输出,需要的正是本章的理论:不可行时去找互相矛盾的约束(Farkas 证书),无界时检查是否漏了预算或仓位约束。 - 哪些风险度量能写成 LP:MAD、最大回撤的某些近似、CVaR(Rockafellar 与 Uryasev 给出的线性化:\(\text{CVaR}_\alpha=\min_\zeta\ \zeta+\frac{1}{(1-\alpha)T}\sum_t(-r_t^Tw-\zeta)^+\),正部同样用辅助变量线性化)。方差不能,它需要二次规划。
- 整数约束让问题变难。持仓只数上限(基数约束)、最小交易单位(整手)、"要么不买、要么至少买 0.5%"这类约束都需要整数变量,问题变成混合整数线性规划(MILP),是 NP 难的。思考题 29-3 的整数间隙说明"先解 LP 再取整"可能明显次优;好在第 26 章的整数性定理告诉我们,网络流结构的 LP 自动有整数解。
- 规模。上例的 LP 有 2080 个变量、1000 多条等式约束,HiGHS 不到一秒。主动收益的约束矩阵是稠密的 \(T\times n\) 块加一个单位阵,现代求解器能利用这种稀疏结构;几千只股票、几年日度数据的 MAD 问题都可以直接求解。
本章小结
线性规划在线性约束下优化线性目标。任何 LP 都可以化为标准型 \(\max c^Tx,\ Ax\le b,\ x\ge0\),再引入松弛变量化为松弛型:基本变量用非基本变量表示,令非基本变量为 0 即得基本解。单纯形法从一个基本可行解出发,每次选目标系数为正的变量入基、用最小比值检验选出基变量、做一次转轴,直到目标系数全部非正;每次转轴 \(O(mn)\),最坏需要指数次(Klee-Minty),但实践中很快。退化可能导致循环,Bland 规则可防止。对偶问题 \(\min b^Ty,\ A^Ty\ge c,\ y\ge0\) 给出原始目标的上界(弱对偶),且最优值相等(强对偶);单纯形法的最终松弛型直接给出对偶最优解——松弛变量系数的相反数,即影子价格。互补松弛刻画了最优性:有富余的约束影子价格为零。辅助线性规划用于判定可行性并获得初始基本可行解。LP 只有三种结局:有限最优、不可行、无界。量化中,L1 跟踪误差、CVaR、线性约束下的 alpha 最大化都是 LP;影子价格给出约束归因和 alpha-跟踪误差前沿的斜率;Farkas 引理就是无套利与状态价格存在性的等价。
| 概念 / 结论 | 公式或要点 |
|---|---|
| 标准型 | \(\max c^Tx\),\(Ax\le b\),\(x\ge0\) |
| 松弛型 | \(z=\nu+\sum_{j\in N}c_jx_j\),\(x_i=b_i-\sum_{j\in N}a_{ij}x_j\)(\(i\in B\)) |
| 基本解 | 非基本变量取 0,基本变量 \(=b_i\) |
| 转轴 | 入基 \(x_e\)(\(c_e>0\)),出基 \(x_l=\arg\min_{a_{ie}>0}b_i/a_{ie}\) |
| 最优性条件 | 目标系数 \(c_j\le0\)(\(\forall j\in N\)) |
| 无界 | 某入基列所有 \(a_{ie}\le0\) |
| 退化 / 循环 | \(b_l=0\) 时目标不变;Bland 规则防循环 |
| 对偶 | \(\min b^Ty\),\(A^Ty\ge c\),\(y\ge0\) |
| 弱对偶 | \(c^T\bar x\le b^T\bar y\) |
| 强对偶 | 最优值相等;\(\bar y_i=-c'_{n+i}\) |
| 互补松弛 | \(\bar x_j(\sum_ia_{ij}\bar y_i-c_j)=0\),\(\bar y_i(b_i-\sum_ja_{ij}\bar x_j)=0\) |
| 辅助 LP | \(\max -x_0\),\(Ax-x_0\mathbf 1\le b\);\(L\) 可行 ⟺ 最优值为 0 |
| 基本定理 | 有限最优 / 不可行 / 无界,三者必居其一 |
| Farkas 引理 | \(Ax\le0,c^Tx>0\) 与 \(A^Ty=c,y\ge0\) 恰有一个有解 |
| 绝对值线性化 | \(\vert e\vert \) → \(u+v\),\(e=u-v\),\(u,v\ge0\) |
练习
基础
- 写出 (29.24)–(29.28) 的 \(n,m,A,b,c\),并给出三个可行解及其目标值。(原书 29.1-1、29.1-2。)
- 把 \(\min\ 2x_1+7x_2+x_3\),s.t. \(x_1-x_3=7\),\(3x_1+x_2\ge24\),\(x_2\ge0\),\(x_3\le0\) 化为标准型。(原书 29.1-4 类型题。提示:\(x_1\) 无符号约束要拆分,\(x_3\le0\) 可令 \(x_3'=-x_3\)。)
- 用 SIMPLEX 手算 \(\max\ 18x_1+12.5x_2\),s.t. \(x_1+x_2\le20\),\(x_1\le12\),\(x_2\le16\),\(x\ge0\),并写出对偶问题和对偶最优解。(原书 29.3-5、29.4-1。可用 29.8.1 的代码核对:最优值 316,\(x=(12,8)\)。)
- 写出最大流 LP(29.47–29.50)的对偶,并解释为什么它对应最小割。(原书 29.4-3。)
- 证明:若在 PIVOT 中目标值不变,则必有 \(b_l=0\)。(原书 29.3-2 相关。)
- 举例说明可行域无界但最优值有限的 LP。(原书 29.1-9。)
进阶
- 证明互补松弛条件(思考题 29-2),并在 29.4.2 节例题上验证。
- 证明 Farkas 引理(思考题 29-4)。提示:考虑 LP \(\max c^Tx\),s.t. \(Ax\le0\) 及其对偶。
- 单变量 LP \(P:\max tx\),s.t. \(rx\le s\),\(x\ge0\)。讨论 \(r,s,t\) 取何值时出现"两者都有有限最优""\(P\) 可行而对偶不可行""对偶可行而 \(P\) 不可行""都不可行"四种情形。(原书 29.5-9。)
- 在 29.8.2 节的跟踪 LP 中加入换手约束:已有持仓 \(w^0\),要求 \(\sum_i|w_i-w_i^0|\le0.2\)。用绝对值线性化写成 LP,求解并报告换手约束的影子价格。
- 用 29.8.2 节的数据实现 CVaR 最小化的 Rockafellar-Uryasev LP:\(\min\ \zeta+\frac{1}{(1-\alpha)T}\sum_tz_t\),s.t. \(z_t\ge-r_t^Tw-\zeta\),\(z_t\ge0\),再加满仓、上限约束,取 \(\alpha=0.95\)。与最小 MAD 组合比较。
- 思考题 29-3:构造一个小的整数规划,使其整数最优值严格小于 LP 松弛的最优值,并计算整数间隙。
原书推荐习题:29.3-5、29.4-2、29.4-3、29.4-5、29.5-5、29.5-9,思考题 29-2(互补松弛)、29-3(整数间隙)、29-4(Farkas 引理)。
原书对照
页码换算:原书页码 = PDF 页码 − 21。
| 本章小节 | 原书章节 | PDF 页码(原书页码) |
|---|---|---|
| 29.1 什么是线性规划 | 第 29 章导言(政治问题、二维例子、算法概述) | p.864–871(843–850) |
| 29.2 标准型与松弛型 | 29.1 Standard and slack forms | p.871–879(850–858) |
| 29.3 把问题写成线性规划 | 29.2 Formulating problems as linear programs | p.880–885(859–864) |
| 29.4 单纯形算法 | 29.3 The simplex algorithm | p.885–900(864–879) |
| 29.5 对偶性 | 29.4 Duality | p.900–907(879–886) |
| 29.6 初始基本可行解 | 29.5 The initial basic feasible solution | p.907–915(886–894) |
| 29.7 思考题与注记 | Problems 29-1 ~ 29-5;Chapter notes | p.915–918(894–897) |