单纯形法¶
单纯形法是解决线性优化问题的高效算法. 虽然它最慢是指数级的, 但它在通常情况下都有很好的性能.
线性优化¶
线性优化问题是指这样的问题:
给定常量矩阵 \(A\), 向量 \(\mathbf{b}\), \(\mathbf{c}\) 与一个非负的变量向量 \(\mathbf{x}\), 求在 \(A\mathbf{x} = \mathbf{b}\) 的情况下, \(\mathbf{c}^T\mathbf{x}\) 的最大值. 取到最大值时的 \(\mathbf{x}\) 称为最优解.
有时, 你会看到约束条件是 \(A\mathbf{x}\leq \mathbf{b}\), 将等于号改成了小于等于. 这种情况下, 只需要为每个不等式增加一个不同的非负变量, 就可以把它转化为等于的情况了. 同样, 如果你的 \(\mathbf{x}\) 并不局限于正数, 而是整个实数, 你可以将每个变量改写为 \(x_1 - x_2\) 的形式, 其中 \(x_1\) 和 \(x_2\) 都是正数. 这两种方法都是通过增加变量的方式把问题转化为一个标准形态, 它们被称作松弛. 下面我们只考虑 \(A\mathbf{x} = \mathbf{b}\) 形式的问题.
我们要求 \(A\) 的行数不大于列数. 对于 \(m\) 行 \(n\) 列的矩阵 \(A\), 其实 \(A\mathbf{x} = \mathbf{b}\) 就是 \(m\) 个等式, 每个不等式含有 \(n\) 个变量. 我们要求等式的数目小于变量数, 否则 \(x\) 就已经完全确定了 (或是根本无解), \(\mathbf{c}^T\mathbf{x}\) 也就随之确定. 我们还要求 \(A\) 满秩, 否则要么 \(\mathbf{x}\) 无解, 要么有一行完全可以消去.
定义 1. 我们把满足 \(A\mathbf{x} \leq \mathbf{b}\) 的 \(\mathbf{x}\) 构成的集合称作一个多面体 (polyhedron). 在线性优化中, 这个多面体被称作可行域 (feasible region).
定义 2. 如果集合 \(S\in \mathbb{R}^n\) 满足 \(\forall x, y\in S, \lambda\in [0, 1], \lambda x + (1-\lambda) y\in S\), 那么 \(S\) 是凸集 (convex set). 换句话说, 如果 \(S\) 中任意两点的连线上的每个点都在 \(S\) 中, 那么 \(S\) 是凸集.
我们把凸集中两点 \(x, y\) 连线上的点称作 \(x, y\) 的凸组合, 也就是上面所提到的 \(\lambda x + (1-\lambda) y\), 其中 \(\lambda \in [0, 1]\).
性质 1. 非空的可行域是凸集.
设 \(\mathbf{x}, \mathbf{y}\) 在可行域中, 依照定义 1, 我们有 \(A\mathbf{x} = \mathbf{b}\), 同时 \(A\mathbf{y} = \mathbf{b}\). 这意味着 \(\forall \lambda\in [0, 1], A(\lambda \mathbf{x} + (1-\lambda) \mathbf{y}) = \mathbf{b}\). 因此, 所有 \(x\), \(y\) 的凸组合 \(\lambda \mathbf{x} + (1-\lambda) \mathbf{y}\) 都在可行域中. 根据定义 2, 可行域是凸集.
定义 3. 凸集中, 无法用两个不同点的凸组合表示的点是顶点.
我们直观意义上理解的"顶点", 确实是不可能在凸集中另外两个点的连线上的. 这个定义只是将"顶点"严格化了.
线性代数告诉我们, \(m\) 行 \(n\) 列的矩阵 \(A\) 所表示的 \(n\) 个变量中, 有 \(m\) 个是线性无关的 (一组基), 它们可以用剩下 \(n-m\) 个变量表示. 让剩下的 \(n-m\) 个变量为 0 所得到的一组解被称为基解. 如果基解在可行域内, 那么就称为基可行解 (basic feasible solution, BFS).
性质 2. BFS 和可行域的顶点一一对应.
先证明 BFS 都在可行域的顶点上. 方便起见, 不妨设前 \(m\) 个变量是线性无关的, 后 \(n-m\) 个均为零. 假设基解 \(\mathbf{x}\) 可以被由可行域中的两点 \(\mathbf{y}, \mathbf{z}\) 凸组合而来, 那么 \(\mathbf{y}\), \(\mathbf{z}\) 的后 \(n-m\) 维必须全部为 0. 这是因为它们的所有变量都非负, 而且 \(\lambda \mathbf{y}_i + (1-\lambda) \mathbf{z}_i = \mathbf{x}_i = 0\). 既然如此, 考虑到 \(A\) 满秩, 剩下的 \(m\) 维只有唯一解. 这说明 \(\mathbf{y}=\mathbf{z}=\mathbf{x}\), 也就是说 \(\mathbf{x}\) 是顶点.
再证明所有的顶点都是基解. 假设顶点 \(\mathbf{x}\) 不是基解, 那么可能存在两种情况: \(\mathbf{x}\) 中等于 0 的维度个数少于 \(n-m\), 或者 \(\mathbf{x}\) 有不多于 \(m\) 个维度非零而且它们线性相关. 第二种情况其实是不可能的.
对于第一种情况, 设顶点 \(\mathbf{x}\) 的 \(m'\geq m\) 个非零维度构成集合 \(I = \lbrace i \mid \mathbf{x}_i > 0 \rbrace\). 那么 \(A\) 的某个列向量 \(A_k\ (k\in I)\) 一定和另外 \(m\) 个列向量 \(A_{k_i}\ (k_i\in I)\) 线性相关 (因为 \(A\) 的秩是 \(m\)). 所以存在一个向量 \(d\), 使得 \(\mathbf{x}_1 = \mathbf{x} + \varepsilon \mathbf{d}\), \(\mathbf{x}_2 = \mathbf{x} - \varepsilon \mathbf{d}\) 都在可行域中 (只需要同时调整那些线性相关的向量就好). 此时, \(\mathbf{x}\) 是它们的凸组合, 所以不是顶点, 矛盾.
性质 3. \(\mathbf{c}^T\mathbf{x}\) 的最大值一定可以在可行域的顶点处取到.
假设这个最大值在可行域内部的某个点 \(\mathbf{y}\) 处取到. 可行域内部的每个点都是所有顶点的凸组合 (想一想, 为什么?). 假设有 \(k\) 个顶点 \(\mathbf{x}_i\ (i = 1, \dots, k)\), 那么
其中 \(\lambda_i\) 的总和为 1. 换言之, \(\mathbf{c}^T\mathbf{y}\) 其实是各个 \(\mathbf{c}^T\mathbf{x}_k\) 的加权平均, 所以这些值中至少有一个不比它小. 因此, 最大值一定可以在某个顶点处取到.
这意味着想要找到最优解, 我们只需要找顶点. 为了找顶点, 我们只需要找到基可行解. 但是, 我们至多有 \(C_n^m\) 个基可行解: 如果 \(A\) 的 \(n\) 列全部线性无关, 从中选出 \(m\) 列作为基解, 将会有 \(C_n^m\) 种选法. 这超过了多项式级复杂度.
我们能否利用某些方法, 快速地选出最优解呢? 这就是单纯形法 (simplex method) 所解决的问题.
单纯形法¶
我们先任意地选择 \(m\) 个线性无关的 \(A\) 的列向量, 构成一组基解. 这组基可以用剩下的 \(n-m\) 个变量表示, 而我们的目标 \(\mathbf{c}^T\mathbf{x}\) 也自然可以用剩下的 \(n-m\) 个变量表示, 得到一个多项式. 不妨记作:
如果这个多项式的系数中, 存在一个 \(\sigma_k\) 是正数, 那么现在的基解不是最优解. 这是因为根据基解的定义, 多项式中的所有变量 \(x_i\) 都等于零; 但如果我们通过更改基解中变量的值, 使得对应的变量 \(x_k\) 不为零, 那么 \(\mathbf{c}^T\mathbf{x}\) 就会变大.
为了进一步靠近最优解, 我们可以将 \(x_k\) 放入基中 (称作替入变量), 同时自然需要将一个变量弹出基 (称作替出变量). 如果有超过一个 \(\sigma_i\) 为正数, 我们会选择最大的那一个替入, 以加快爬升向最大值的速度. 同理, 我们会选择"弹出之后, 使得新加入的变量最大"的变量替出基, 这也是为了加快爬升的速度. 我们重复进行这样的替换, 直到所有系数都为 0 或负数为止.
接下来我们只需要解决两个问题:
- 计算弹出一个基变量, 对当前选中的替入变量产生的影响;
- 计算加入一个基变量后, \(A\) 的改变与 \(\mathbf{c}^T\mathbf{x}\) 的新系数.
对于第一个问题, 在求出基解后, 矩阵 \(A\) 应当经过了初等行变换, 使得基解所在的列向量构成一个单位矩阵. 不妨假设我们要替入 \(x_k\), 它对应的列向量是 \(A_{\cdot k}\). 这意味着, 如果第 \(i\) 个基解 \(x_{B_i}\) 降低 1, \(x_k\) 就应当升高 \(1 / A_{ik}\), 以保证那一行的等式仍然成立. 所以, 我们只需要求出:
所得到的 \(i\) 就指示着应当被替换出去的基变量 \(x_{B_i}\). 当然, 我们需要保证替换完成之后所有的基变量都是正的.
对于第二个问题, 依然假设替入变量是 \(x_k\), 对应的列向量是 \(A_{\cdot k}\), 而它要变成第 \(j\) 个基解.
首先, 我们将这个向量除以 \(A_{jk}\), 也就是使得 \(A_{jk}\) 恰好为 1. 这时, 等式 \(A\mathbf{x} = \mathbf{b}\) 仍然成立, 这是因为 \(x_k\) 此时仍然是非基变量, 等于零.
接下来, 对每个 \(i\neq j\), 将第 \(i\) 个基解 \(x_{B_i}\) 所在的列乘以 \(A_{ik}\), 然后令 \(x_k\) 所在的列减去它, 就可以使那一列的第 \(i\) 个元素变为零. 当然, 等式的另一端也要如此处理.
那么 \(\mathbf{c}\) 怎么办呢? 实际上, 我们在执行单纯形法的时候计算 \(A\mathbf{x}\leq \mathbf{b}\) 的约束条件时, 采用的是如下的扩充矩阵:
通过引入单位矩阵 \(I\), 我们将不等式松弛为等式, 而且可以免去最开始求出基解的高斯消元. 通过在合适的位置加入不等式的右端 \(\mathbf{b}\) 和目标函数 \(\mathbf{c}\), 我们在更新 \(A\) 的系数时, 也一并更新了它们.
然而，这种解法要求 \(b>0\), 否则最开始的基解就不是可行解. 这是因为此时基解就恰好等于 \(\mathbf{b}\). 为了解决这个问题, 我们可以参考对偶单纯形法.
这就是单纯形法的全貌. 具体的实现可以参看 OI Wiki.