线性规划单纯形法代码实现解析
本文仅针对代码实现思路，并不完整讲解单纯形算法原理，建议先完成课本或资料的原理学习，再阅读本文。
1. 算法课本的实现
本节代码完全按课本方式实现，竞赛选手可略过此部分，查看第2节。
1.1 流程梳理
首先，我们先回顾单纯形基本流程：
- 确定初始基本可行解
- 检查检验数 \(\lambda\) 和 \(B^{-1}A\) 的 \(\alpha\) 矩阵做判断：
- 是最优解或无最优解则结束；
- 否则，做基变换，合理选择一个非基变量换入，一个基变量换出，更新变换后的各值
- 重复2
1.2 参数定义
先看看课本的单纯形法有哪些量：
- \(m\) ：约束条件的行数
- \(n\) ：变量的个数
- \(A\) 矩阵： \(m\times n\) 的矩阵，约束条件的系数矩阵，通常情况 $ m < n $ 且 \(A\) 的秩为 \(m\) 。如果完整构建，则该矩阵应该也存了松弛变量对应的单位矩阵那部分
- \(b\) ：约束条件右侧的值向量，最小标准形应该不小于0
- \(c\) ：目标函数的系数向量，最小标准形就是求 \(\min{c^{T}x}\)
- \(z\) ：目标函数，不过表达上经过了变换，在算法中每次由上一次的目标函数值来计算
- \(B\) 基： \(m\times m\) 的矩阵，为 \(A\) 中选出的 \(m\) 列线性无关的一组基，如果是完美的松弛形，则恰好可以选由松弛变量系数构成的“最右边”那个 \(m\times m\) 的单位矩阵
- \(\alpha=B^{-1}A\) 矩阵：由基 \(B\) 和约束矩阵 \(A\) 得到的校验参数
- \(\beta=B^{-1}b\) ：其实就是枚举的每一组基本可行解
我们来个课本例子：
\[ \begin{align} \min ~~ & {z=-12x_{1}-15x_{2}} \\ s.t. ~~ & 0.25x_{1} +0.5x_{2}& + x_{3} &&& =120 \\ & 0.5x_{1} +0.5x_{2}&& + x_{4} && =150 \\ & 0.25x_{1}&&& + x_{5} &=50 \\ & x_{i}\geq 0, i=1,2,\dots, 5 \end{align} \]
用表来展现：
| -12 | -15 | 0 | 0 | 0 | c/b |
|---|---|---|---|---|---|
| \(x_1\) | \(x_2\) | \(x_3\) | \(x_4\) | \(x_5\) | |
| 0.25 | 0.50 | 1 | 0 | 0 | 120 |
| 0.50 | 0.50 | 0 | 1 | 0 | 150 |
| 0.25 | 0 | 0 | 0 | 1 | 50 |
\(A\) 矩阵就是完整的系数矩阵：
| 0.25 | 0.50 | 1 | 0 | 0 |
| 0.50 | 0.50 | 0 | 1 | 0 |
| 0.25 | 0 | 0 | 0 | 1 |
第一组可行基\(B\)我们就选右侧的单位矩阵：
| 1 | 0 | 0 |
| 0 | 1 | 0 |
| 0 | 0 | 1 |
既然选的这组基对应的变量序号是 \(3,4,5\) 列，那么基变量 $x_{B}=x_3,x_4,x_5$ ，非基变量 $x_{N}=x_{1},x_{2}$ 。
相应的，基变量对应的目标函数的系数 $c_{B}=,0,0$ ，非基变量对应的目标函数的系数 $c_{N}=,-15$ ，发现了吗，松弛变量因为是我们构建标准形额外加的，所以它们序号对应的目标函数的系数都是0，也即最初的基变量对应的目标函数的系数都是0 。
记住以上内容（或文章的位置），我们接下来的代码实现会参照这些内容。
1.3 搭个框架
1 | struct Simplex |
这里，我们用a[][]、b[]、c[] 保存最初的 \(A\) 矩阵、 \(b\) 、 \(c\) 向量。
1.4 实现
1.4.1 循环枚举可行基
算法中“重复2”的过程体现在Solve()的while(true)里，在这里检查是否结束。 如果需要基变量置换，则在Solve()中调用Pivot()执行置换过程。
Solve()的开头，参照课本最优性检验的定义， \(k\) 是换入变量序号， \(l\) 是换出变量在基里的序号，先把这俩序号定义出来。
1 | double Solve() |
如果 \(\lambda \geq 0\) 则是最优解直接返回，否则如果有 \(\lambda_k < 0\) ，则 \(k\) 就是换入变量的序号。
代码用一个单行for循环找到 \(k\) ，利用了for循环里的判断条件，注意这个for的结尾直接跟分号“;”。如果循环结束k==n的话那就是没找到小于0的 \(\lambda\) ，这C语言基础知识了，这种情况已经最优解，属于“情况1”，返回当前的目标函数值z。
1 | for(k = 0; k < n && c[k] >= 0; k ++); |
困惑来了：不是 \(\lambda\) 吗，这里怎么用c[k]>=0做条件？我们回顾课本检验数的公式： \(\lambda^{T}=c^{T}-c_{B}^{T}B^{-1}A\) ，翻一下上文的 \(c_{B}\) 和 \(c_{N}\) ，当松弛变量为初始基的时候，对应的 \(c_{B}\) 都是0呢，所以“恰好” \(c\) 就是初始的 \(\lambda\) 。由于单纯形法的过程中都是由“上一步”的 \(\alpha,\beta,\lambda,z\) 推“下一步”的，所以初始的 \(c\) 可以不用保留，我们直接用c[]数组表示接下来所有的 \(\lambda\)。
接下来就处理有那么个 \(\lambda_{k}<0\) 的情况，看 \(\alpha\) 的第 \(k\) 列是否有大于0的。
1 | double mn = inf; |
又有心细的同学问了，不是 \(\alpha\) 吗，怎么直接用a[][]矩阵了？回忆上文 \(\alpha\) 是什么？是 \(B^{-1}A\) ，而初始基 \(B\) 是单位矩阵，所以初始的 \(\alpha\) 恰好就是 \(A\) 矩阵，我们又可以直接用 a[][]表示后续的 \(\alpha\) 矩阵了。
如果临时变量mn==inf，意味 \(\alpha\) 的第 \(k\) 列没有遇到大于0的，根据算法，这种情况无最优解，属于“情况2”，返回一个信号，这里我们从简，返回个无穷大好了。
for循环里面就是在顺便处理“情况3”了，即确定换出变量的序号 \(l\) ，根据算法，是 \(\alpha\) 这第 \(k\) 列大于0的里面，对应的 \(\beta_{i}/\alpha_{i,k}\) 最小的那个 \(i\) 作为 \(l\) 。
如果“情况1”和“情况2”都没有遇到，即函数没有return，就要根据找到的 \(k\) 和 \(l\) 做基变换了，我们调用Pivot(k,l)，完整的Solve()代码为：
1 | double Solve() |
1.4.2 基变换
接下来看void Pivot(int k, int l)的实现。
回顾课本公式：
- \(\alpha_{lj}^{\prime}=\alpha_{lj}/\alpha_{lk}, 1\leq j \leq n\)
- \(\alpha_{ij}^{\prime}=\alpha_{ij} - \alpha_{ik}\alpha_{lj}/\alpha_{lk}, 1 \leq i \leq m 且i\neq l, 1\leq j \leq n\)
- \(\beta_{l}^{\prime}={\beta_l}/\alpha_{lk}\)
- \(\beta_{i}^{\prime}=\beta_{i} - \alpha_{ik}\beta_{l}/\alpha_{lk}, 1\leq i \leq m 且 i \neq l\)
- \(\lambda_{j}^{\prime}=\lambda_{j}-\lambda_{k}\alpha_{lj}/\alpha_{lk}, 1\leq j \leq n\)
- \(z_{0}^{\prime}=z_{0}+\lambda_{k}\beta_{l}/\alpha_{lk}\)
我们就依据公式把基变换之后所有的数都更新出来。
先看l这一行，这里我们一开始没有更新 \(\alpha[l][k]\)（即a[l][k]），因为它后面还要被用到，最后再更新它。
1 | b[l] /= a[l][k]; |
接着处理除了第l行之外的每一行，a[i][k]==0时候可以跳过不用计算，看公式就知道为0时对结果没有影响。
1 | for(int i = 0; i < m; i ++) |
这里更新a[i][j]时的代码为a[i][j] -= a[i][k] * a[l][j];，公式是 \(\alpha_{ij}^{\prime}=\alpha_{ij} - \alpha_{ik}\alpha_{lj}/\alpha_{lk}\) ，有同学可能要问“为什么少除了 a[l][k] 呢”？因为前面的代码已经更新了 a[l][j]的值，到这行代码时的 a[l][j] 已经是 \(\alpha_{lj}/\alpha_{lk}\) 了。
这里我们对第k列做了单独处理（置0），一方面计算过程中要用到a[i][k]的值，不能过早覆盖它，另一方面其实按公式算也是0，浮点数运算直接置0也能优化点精度。
接着更新目标函数 \(z\) 和检验数 \(\lambda\) ，前面已经讲解了我们直接用c[]数组保存后续更新的 \(\lambda\)。
1 | z += c[k] * b[l]; |
c[k]单独处理理由和前面一样的，过程中要用它不能过早覆盖。
最后别忘了前面的a[l][k]也要更新，按公式就是除以它自己，应该为1。Pivot()完整代码如下：
1 | void Pivot(int k, int l) |
1.4.3 初始化
代码中要用到行数m和列数n，当然要初始化好，可以定义个内置函数
1 | void Init(int n_, int m_) |
重点在于所有a[][]、b[]、c[]的值要根据输入初始化，不要在处理多组数据时残留之前的数据，初始化为0比较安全。vector容器增强代码的可伸缩性。
当然，习惯用C语言形式的数组也可。
1 | const int maxn = 1e3 + 10; |
特别注意的是，当存入了松弛变量，实质上变量数量增加了，传入Init()的n要对应哦。
1.4.4 使用
以P3980 [NOI2008] 志愿者招募 这道题为例：
n天，第i天需要bi个人- 有
m类志愿者，第j类可以从第si天工作到第ti天，每人ci元- 第一行输入
n m，第二行n个数表示每天需要人数，接下来m行每行 3 个数表示一类志愿者的起止时间和花费- 求满足要求的最省钱总费用
直接建线性规划模型，设x是每类志愿者个数，约束条件每一行是每一天的要求，\(A\)矩阵第i行为01向量，\(A_{ij}\) 表示第j类志愿者第i天是否可以服务。
\[ \begin{align} \min ~~ &{z=cx} \\ s.t. ~~ &\sum A_{ij}x_{j} \geq b \\ &x\geq 0 \end{align} \]
这全是 \(\geq\) ，不便于我们建松弛变量，转成对偶问题来做：
\[ \begin{align} \max ~~ & {z=by} \\ s.t. ~~ & \sum A_{ji}y_{i} \leq c \\ & y\geq 0 \end{align} \]
这还不够，不是最小标准型，把目标函数系数取个负号变成
\[ \min{z=-by} \]
求解完成后记得把负号乘回来。
接下来实现构造数据的代码：
1 | Simplex spx; |
给一份完整代码：
1 |
|
“啊！怎么MLE了，这不是正解啊？！”
是的，完全按课本来，我们遇到了一个问题，加了松弛变量之后，数据的两个维度都达到了这道题m 的 \(10^4\) 级，a[][]矩阵达到 \(10^8\) ，虽然过了部分数据，但大数据 Memory Limit Exceed 了。
看来这个写法只能解决规模没那么大的问题，那么应该怎么办呢？请看下一节。
2. 改进的实现
2.1 松弛变量与置换策略
对于刚刚的例子，n的范围在 \(10^3\) ，m的范围在 \(10^4\) ，如果只是 \(n\times m\) ，还不至于MLE，但加上松弛变量，a[][]就要开到 \((n+m)\times m\) 。
我们知道松弛变量系数构成的矩阵是个单位矩阵，能不能不显式地保存它，却仍利用它为初始的可行基 \(B\) 进行计算呢？
在构造数据的时候，先不保存松弛变量，注释掉的即我们调整的地方：
1 | // spx.Init(n + m, m); --> |
然后需要对Pivot()进行改变：
1 | void Pivot(int k, int l) |
我们改变了三个地方，对应着特殊的k 这一列。课本上利用 \(H\) 矩阵把 \(\alpha_{lk}\) 变为1，k这列的其他 \(\alpha_{ik}\) 消元为0，对应的 \(\lambda_{k}\) 也变为了0。
而改变之后，a[][]数组里没有再存储松弛变量对应的单位矩阵，基变换的过程中也没有了那些列的信息，第k列相应改变了基变换时更新的方式，a[l][k]变为1/a[l][k]，其它的a[i][k]变为-a[i][k]/a[l][k]。
“不对呀，这跟课本的公式就不一样了啊，为什么按课本的方法要在a[][]里存松弛变量，而这个不用存？不用存就罢了，a[][k]这第k列为什么这样更新？”
我们需要先理清课本方法做了什么事：
把a[][]的第k列的a[l][k]变为1，其它a[i][k]变为0，是 \(H\) 矩阵的作用，即让可行基 \(B\) 乘特意设计的 \(H\) 矩阵完成基变换，将原始 \(A\) 矩阵的第 \(k\) 列换入，基 \(B\) 的第 \(l\) 列换出，导致的 \(\alpha\) 矩阵的变化。我们需要存储松弛变量对应的那些列，才能完整体现换入和换出的所有信息。
以这样一组数据为例：
1 | 10 10 |
我们以课本方法进行一次换入换出，打印一个中间步骤的l、k、a[][]( \(\alpha\) )和c( \(\lambda\) )：
1 | l: 7 k: 9 |
经过公式计算，第12列的c[]是-4，这里加上了中括号[]方便找到。这一列在a[][]中一开始存储的是松弛变量的系数。
有一点不要混淆：我们只是因为初始的a[][]恰好等于 \(\alpha\) 矩阵，才用a[][]数组继续保存后续的 \(\alpha\) ，但 \(A\) 矩阵和 \(\alpha\) 矩阵并不等价， \(A\) 是不会变的，改变的是完成变量置换这个任务的 \(H\) ，且 \(\alpha=H^{-1}B^{-1}A\) 。
如果a[][]中不保存松弛变量系数，则这些列便不再存在，在更新 \(\alpha\) 的过程中，这一次c[]的第12列为负的信息就无从知晓，本应当在某一次迭代作为换入变量的它却遗失了，导致后续计算错误。
知道这个情况之后，我们再考虑改进的方法，不保存松弛变量，那么a[][]的第k列新的更新方式是什么原理呢？我们打印一份第一步的结果，对比课本方法带松弛变量系数列的结果，和改进方法不带松弛变量列的结果：
课本方法：
1 | l: 9 k: 0 |
改进方法：
1 | l: 9 k: 0 |
更新的是换入变量k=0列（数组下标从0开始），换出变量l=9列（最后一列）。
我们看到，课本方法最后一列的 \(\lambda\) 变为了 8（用中括号[]标记了），换入变量的k=0的 \(\lambda\) 变为了0 。
于是发现，改进方法并不是基于 \(H\) 矩阵更新第k=0列的值，而是把k=0列变成了更新后的l=9那一列的 \(\alpha\) 值。
这样做，就保留了课本方法里松弛变量对应的那些列的更新信息，不会缺失后期需要变量置换的列的信息。
“那原本k=0那一列就丢掉了？”
是的，暂时“丢掉了”，它被“换入”了，一个已经在基中的变量，不可能再次被选中进行换入，所以暂时丢掉它的信息，不影响我们后续寻找用于换入的变量。
这样每次都用换出变量的列序号（\(l\)）对应的 \(\alpha\) 值替换换入变量对应的列（\(k\)）的 \(\alpha\) 值，就使得课本方法里存储的松弛变量那些列序号被换入后的信息，存在了一个“虚空”里，它会基于数学推导的结果被换入换出，而不影响计算结果。
2.2 检验数的选取
肯定会有多个 \(\lambda\) 分量小于0的情况，无论选哪个都肯定对，这次没换入，下次也就换入了。
但是如果更深入了解单纯形与凸优化的话，会发现，选取最小的负\(\lambda_k\)作为换入变量，会让算法收敛更快，起到提速作用。
1 | double Solve() |
2.3 题外话：浮点数计算的好习惯
浮点数在计算过程中有精度损失，一个本应为0的数，就可能在0的左右浮动，如果我们直接判断它是否大于等于0，而恰好浮动在很小的-0.00...01，就得到了错误的判断。
解决方法是定义一个全局的精度控制\(eps=10^{-8}\)，这个大小只是经验之谈，\(10^{-7}\)、\(10^{-6}\)也不是不行。
这时我们判断一个数是否大于0，就将if(x > 0) 改为 if(x > eps)，如果判断是否大于等于0，就将if(x >= 0) 改为 if(x > -eps)。
两数比大小也类似，if(x > y) 改为 if(x - y > eps)。
2.4 以例题为例的完整代码
该代码可做模板保存使用。
1 |
|
3. 输出解
如果要输出求解 \(x\) 向量的值，我们需要记录换入换出的情况，增加两个数组IdA[], IdB[] 分别记录原始变量的序号id、基变量的序号id，在Pivot()函数增加换入换出后的序号信息：
1 | void Pivot(int k, int l) |
如何初始化IdA与IdB，最终如何输出解，就留给大家尝试吧。
本文不涉及单纯形的原理和数学证明，仅解析代码实现。
如有原理性错误欢迎指正。
主要参考：
- 屈婉玲，刘田，张立昂，王捍贫，《算法设计与分析》，清华大学出版社，2014.
- 网上的单纯形模板