单纯形算法详细解析

线性规划（Linear Programming，LP）是很是经典的算法之一，而解决该问题的最经常使用方法是单纯形法。本博文致力于用最简单、最详细的语言一步步解释单纯形算法的过程并加以详细的解释。

中学课程里，咱们都简单地接触过线性规划，那时候通常都是分析每一个约束，在二维平面上画出直线，获得可行域，而后以固定斜率做出目标函数直线，在可行域内移动直线，在y轴上的截距就是最优解。而每每最优解的地方是经过（凸）可行域的顶点。就像下面这个例子：
\[ \begin{equation} \begin{split} \max& \quad x_3+x_4 \\ s.t.& \quad -x_3+2x_4\leq 2 \\ &\quad 3x_3-2x_4\leq 6 \\ &\quad x_3, x_4\geq 0 \end{split} \end{equation} \]

蓝色区域是可行域，红色直线是固定斜率的，当上移到（4，3）点时目标函数取最大值。

而咱们后面将证实，最优解必定是可行域（凸超几何体）的顶点之一。那么，咱们先假定这成立，就使用”改进-中止“（improve-stopping）的套路，即给定可行域的一个顶点，求值，沿一条边到达下一个顶点，看是否能改善解，直到达到中止要求。

这里就要问几个问题了：为何最优解必定在顶点处？怎么获得顶点？怎么实现沿着一条边到下一个顶点？何时中止？接下来，咱们将一一解答这些问题，当解答完这些问题，单纯形法也就显而易见了。

一、凸多边形最优解在顶点

考虑最小化问题，目标函数\(\mathbf{c}^T \mathbf{x}\)，有一个在可行域内部的最优解\(\mathbf{x}^{(0)}\)，由于凸多边形内部任一点均可以表示成顶点的线性组合，即对于顶点\(\mathbf{x}^{(k)}, k=1,2,\cdots, n\)，有
\[ \mathbf{x}^{(0)}=\sum_{k=1}^{n}\lambda_k\mathbf{x}^{(k)}, \]

\[ \sum_{k=1}^n \lambda_k=1 \]

假设\(\mathbf{x}^{(i)}\)是全部顶点中使得\(\mathbf{c}^T \mathbf{x}\)最小的顶点，那么有
\[ \begin{equation} \begin{split} \mathbf{c}^T\mathbf{x}^{(0)} &= \sum_{k=1}^{n}\lambda_k\mathbf{c}^T\mathbf{x}^{(k)} \\ &\geq \sum_{k=1}^{n}\lambda_k\mathbf{c}^T\mathbf{x}^{(i)} \\ &= \mathbf{c}^T\mathbf{x}^{(i)} \end{split} \end{equation} \]
所以，总有一个顶点，他的目标函数值不比内部点差。

二、怎样获得一个凸多边形的顶点？

下面要说明的就是这样一个定理：对于可行域方程组\(\mathbf{Ax=b}\)，该方程肯定的凸多边形的任意一个顶点对应\(\mathbf{A}\)的一组基。

2.1 顶点对应一组基

上面这个例子是松弛形式的约束，原来的变量有三个，加上后面4个后变成等式，造成的可行域如上图所示。咱们取一个顶点（0，2，0）分析，代入约束中，能够算出一个完整解\(x=(0，2，0，2，2，3，0)\)。取出矩阵\(\mathbf{A}\)中对应的\(x\)不为0的列，即图中标蓝的几列（用\(\mathbf{a}_2\)，\(\mathbf{a}_4\)，\(\mathbf{a}_5\)，\(\mathbf{a}_6\)表示），那么这几列就是线线性无关的，考虑\(m<n\)（约束数目小于松弛后的变量个数），那么自由解有\(n-m\)维，所以挑出来的列必有\(m\)个，构成一组基。下面主要说明他们为何线性无关。假设他们线性相关，即存在一组\(\lambda\neq\mathbf{0}\)，使得\(\lambda_2\mathbf{a}_2+\lambda_4\mathbf{a}_4+\lambda_5\mathbf{a}_5+\lambda_6\mathbf{a}_6=0\)，也就是说，能够构造\(\lambda=[0, \lambda_2, 0, \lambda_4, \lambda_5, \lambda_6, 0]\)，使得\(\mathbf{A}\lambda=0\)。那么还能够再构造两个异于\(x\)的解：\(x'=x+\theta\lambda\)和\(x''=x-\theta\lambda\)。他们都知足\(\mathbf{Ax=b}\)。而且能够经过控制\(\theta\)取很小的正值，使得这两个解知足都大于0的约束。由此这两个解都在凸多面体内，而且有\(x=\frac{1}{2}(x'+x'')\)，可是这是有问题的，由于一个凸多面体的顶点是不能被内部点线性表示的，所以这几列是构成一组基的。

在这里，咱们还能够对每一组解，都将 \(\mathbf{A}\)的列从新排列一下，将解向量也排列一下，写成分块矩阵的形式，那么就会有 \(\mathbf{x}_{\mathbf{B}}=\mathbf{B}^{-1}\mathbf{b}\)和 \(\mathbf{c^Tx=c_B^TB^{-1}b}\)。这是两个颇有用的式子，在后面单纯形算法的理解上颇有帮助，这里先记下。

2.2 一组基对应一个基可行解（顶点）

有了上面的知识，给定一组基\(\mathbf{B}\)，咱们直接构造\(\mathbf{x}=[\mathbf{B^{-1}b}, \mathbf{0}]^T\)，只要说明他不能被两个异于他的内部点线性表示便可。假设有两个内部点\(\mathbf{x}_1, \mathbf{x}_2\)，知足\(\mathbf{x}=\lambda_1\mathbf{x}_1+\lambda_2\mathbf{x}_2\)，因为\(\mathbf{x_N}=0\)，而且这些解的元素都大于等于0，所以\(\mathbf{x_1}_{\mathbf{N}}=\mathbf{x_2}_{\mathbf{N}}=0\)。而又由于\(\mathbf{Ax_1=Ax_2=b}\)，所以\(\mathbf{x_1}_{\mathbf{B}}=\mathbf{x_2}_{\mathbf{B}}=\mathbf{B^{-1}b}\)。即这两个解和\(\mathbf{x}\)相同，所以\(\mathbf{x}\)是顶点。

3 如何从一个顶点沿着边到另外一个顶点？

这里是要研究怎么改善一个解，咱们须要知道怎么从一个顶点出发沿着边找到另外一个顶点。前面咱们知道了一个顶点对应一组基，并且一个矩阵的基有多个，那么是否能够经过基的变换从而使得顶点变换呢？先来看一个例子。

图中红色点对应一个彻底解 \(\mathbf{x}=[2,0,0,2,0,3,6]\)，对应的基是 \(\mathbf{B}=\{\mathbf{a}_1，\mathbf{a}_4，\mathbf{a}_6，\mathbf{a}_7\}\)，考虑向量 \(\mathbf{a}_3\)，即绿色那列，他能够表示成：
\[ \mathbf{a}_3=0\mathbf{a}_1+1\mathbf{a}_4+1\mathbf{a}_6+1\mathbf{a}_7 \]
将式子补全，就会有
\[ 0\mathbf{a}_1+0\mathbf{a}_2-1\mathbf{a}_3+1\mathbf{a}_4++0\mathbf{a}_5+1\mathbf{a}_6+1\mathbf{a}_7=0 \]
把系数写出来： \(\lambda=[0,0,-1,1,0,1,1]\)，他就是对应上图中的绿色向量（相反方向）。所以，只有沿着这个方向走适合的步长 \(\theta\)，就能到达下一个顶点。即新的顶点和旧的顶点关系为：
\[ \mathbf{x'}=\mathbf{x}-\theta\mathbf{\lambda}(\theta>0) \]
那么多大的 \(\theta\)合适呢？咱们很容易知道 \(\mathbf{x'}\)可以知足 \(\mathbf{Ax=b}\)，由于 \(\mathbf{A\lambda=0}\)，如今要保证的就是 \(\mathbf{x'}\)的各个份量大于等于0。对于 \(\lambda_i\leq0\)的项，相减后大于0，没问题。可是对于 \(\lambda_i>0\)的项，就要当心了，为了保证相减后仍然大于等于0，咱们设置
\[ \theta=\min\limits_{\mathbf{a}_i\in \mathbf{B},\lambda_i>0}\frac{x_i}{\lambda_i} \]
就能保证 \(\mathbf{x'}\geq0\)。在上面的例子中， \(\theta=2\)，新解是 \(\mathbf{x'}=[2,0,2,0,0,1,4]\)，对应的基是 \(\mathbf{B'}=\{\mathbf{a}_1，\mathbf{a}_3，\mathbf{a}_6，\mathbf{a}_7\}\)，到这里，一切看上去很完美，咱们找到了运动到下一个顶点的方法，也就是先找一个非基向量，将他写成用基向量表示的形式，提出系数，肯定步长，获得新解。可是还有一个小问题，咱们看到实际上 \(\mathbf{B'}\)和 \(\mathbf{B}\)差了一个向量，至关于把 \(\mathbf{a_4}\)换出去，把 \(\mathbf{a}_3\)换进来。咱们称这个过程为换基，后面算法实现部分叫pivot。那么，怎么保证换了个向量以后，仍然是基呢？证实一下：

证实：\(\mathbf{B'}=\mathbf{B}-\{\mathbf{a}_l\}+\{\mathbf{a}_e\}\)仍然是基。（l表示leave，e表示enter）

假设\(\mathbf{B'}\)线性相关，那么存在\(<d_1,d_2,\ldots,d_{l-1},d_{l+1},\ldots,d_m, d_e>\neq0\)，使得\(\sum_{k}d_k\mathbf{a}_k=0\)。而\(\mathbf{a}_e=\sum_{i=1}^m \lambda_i\mathbf{a_i}\)，代入得：
\[ (d_1+d_e\lambda_1)\mathbf{a_1}+\ldots+(d_e\lambda_l)\mathbf{a}_l+\ldots+(d_m+d_e\lambda_m)\mathbf{a}_m=0 \]
这里是证实的关键之处：咱们在设置\(\theta\)时的作法，假如最终选出来的使得\(\frac{x_i}{\lambda_i}\)最小的那一项下标为\(p\)，那么获得的新解的第\(p\)项必然为0，至关于把\(\mathbf{a}_p\)换了出去。在上面这个例子重\(p=4\)。而由于咱们只考虑\(\lambda_i>0\)的基向量，所以被换出去的基\(\mathbf{a}_l\)对应的\(\lambda_l>0\)，所以上式中有\(d_1=d_2=\ldots=d_m=d_e=0\)，和原假设矛盾，所以\(\mathbf{B'}\)也是线性无关的。

截止到目前，咱们回答了三个问题，基本上单纯形算法呼之欲出了，单纯形算法就是经过反复的基变换（经过向量的进出）来找顶点，从而找到达到最优值的顶点。可是仍是有些细节须要琢磨，好比，怎么选入基向量？改善的过程何时中止？

4 入基向量的选择及中止准则

以最小化问题为标准，咱们的最终目标是最小化\(\mathbf{c^Tx}\)，所以一个很天然的贪心想法是每步的改善都尽量地大，所以能够计算一下更新的目标函数值和原来的差值，取使得变化大的顶点继续下一步迭代。那么这个差值怎么可以向量化地计算呢？只有向量化地计算，才能避免一个一个地计算比较，提升效率。

假设\(\mathbf{B}=\{\mathbf{a}_1, \mathbf{a}_2,\ldots, \mathbf{a}_m\}\)，那么对于任何一个非基向量\(\mathbf{a}_e\)，都有\(\mathbf{a}_e=\lambda_1\mathbf{a}_1+\ldots+\lambda_m\mathbf{a}_m\)。将\(\lambda\)写完整：\(\lambda=[\lambda_1,\lambda_2,\ldots,\lambda_m,0,0,\ldots,-1,\ldots,0]^T\)，差值\(\mathbf{c^Tx'-c^Tx=c^T(-\theta\lambda)=\theta(c_e-\sum_{a_i\in B}\lambda_ic_i)}\)，所以咱们要选使得这个值的绝对值最大的\(\mathbf{a}_e\)。那么何时表示找到最优值应该中止呢？很明显，就是对于全部\(\mathbf{a}_e\)，这个差值都大于等于0，即目标函数再也不减少。所以，每次迭代，先计算差值，若是存在小于0的，就选一个使得差值绝对值最大的做为入基向量。

嗯，接下来就是要考虑向量化操做。首先咱们看一下\(\lambda\)能不能向量化表示：因为\(\mathbf{B}\lambda=\mathbf{a}_e\)（\(\lambda\)只取基系数部分），所以\(\lambda=\mathbf{B^{-1}}\mathbf{a}_e\)，若是对整个矩阵\(\mathbf{A}\)左乘\(\mathbf{B^{-1}}\)，这就颇有意思了，全部的非基列将变成该非基向量用基向量表示的系数，而全部的基列将变成\(e_k\)，即合起来成为一个单位阵的形式。这是很关键的一步，在单纯形算法实现中也涉及到，先记下。进一步，咱们取\(\mathbf{c}\)的基部分\(\mathbf{c_B}\)，这样\(\mathbf{c_B^TB^{-1}A}\)就等于了上式中的\(\sum_{\mathbf{a}_i\in \mathbf{B}}\lambda_i\mathbf{c}_i\)的向量化形式（对全部的非基向量一同操做）。而后再加上一部分，变成\(\mathbf{c^T-c_B^TB^{-1}A}\)，这就是最终的形式，称之为检验数\(\bar{\mathbf{c}}\)。很容易验证，基向量对应的检验数都是0，咱们的目标就是经过迭代，使得\(\bar{\mathbf{c}}\geq0\)，这时对于任何一个可行解\(\mathbf{y}\)（\(\mathbf{Ay=b,y\geq0}\)），都有\(\mathbf{c^Ty}\geq\mathbf{c_B^TB^{-1}Ay=c_B^TB^{-1}b=c_B^Tx_B=c^Tx}\)，即\(\mathbf{x}\)就是最优的。

5 单纯性算法核心：单纯形表

终于，一系列的讲解结束，单纯性算法也就瓜熟蒂落了。要将上面的一堆东西整理成一个简短高效的可行算法并不简单。先贴上算法伪代码：

再献上一个很是漂亮的单纯形表：

如今，咱们来对算法进行分析，将算法的每一个操做和咱们上面的介绍对应起来。

• SIMPLEX算法一开始调用INITIALIZESIMPLEX找到一个初始基可行解，这在某些状况下很简单，当\(\mathbf{b}\geq0\)时，他就是一个初始基可行解，不然，可能要用到两阶段法、大M法等求，这不是重点。
• WHILE循环内，只要找到一个\(c_e<0\)，就继续迭代。第10到16行就是经过设定\(\theta\)找到出基向量。
• 最关键最有意思的就是PIVOT算法，他巧妙地将咱们介绍的繁杂操做使用一个简单的高斯行变换就实现了。而这个算法的载体就是单纯形表，如上图，左上角是目标函数值相反数\(-z\)，第一行是检验数\(\bar{\mathbf{c}}\)，左下角是基对应的部分解（其余部分是0，不用写出），右下角是矩阵\(\mathbf{A}\)。他始终被分块成两部分，基向量部分始终以单位阵的形式存在。注意左边的部分解每一个份量都是严格对应着一个基向量，即他们是有顺序的。
• 一行一行地看PIVOT算法。输入参数告诉咱们下标为\(l\)的向量被换出，所以找到他对应的那行，暂称为第\(l\)行，这一行对应的基的下标要被换成\(e\)，那么为何更新后对应的解是\(\frac{b_l}{a_{le}}\)呢，要注意，其实这个值就是\(\theta\)，\(b_l\)就是旧的\(x_i\)，\(a_{le}\)就是\(\lambda_i\)（上面已经解释了乘上\(\mathbf{B^{-1}}\)后每一列都是系数）。那么为何更新后是\(\theta\)呢？咱们回到式子\(\mathbf{x'=x-\theta\lambda}\)，因为\(b_l\)如今对应的新向量不是\(\mathbf{x}\)对应的基向量，所以\(\mathbf{x}\)在该位置的值是0，而咱们知道\(\lambda\)在入基向量对应的位置的值是-1，所以\(0-(-1)\theta=\theta\)。
• 第3到4行，将第\(l\)行除以\(a_{le}\)，目的就是将\(a_{le}\)变成1，由于要始终保持基是以单位阵的形式存在。
• 第8行，就是在执行\(\mathbf{x'=x-\theta\lambda}\)的操做，获得新解。
• 第10行，高斯行变换，你会发现这样操做完后，入基列就变成和刚才出基列同样，高斯行变换保证了矩阵的性质。
• 第14行，咱们知道\(-z=0-\mathbf{c_B^TB^{-1}b}\)，因为旧有的基对应的\(c\)都是0，而只有新换入的向量对应的\(c_e\)不为0，具体写一下，减掉的那部分就只有\(c_e\)和他对应的解\(b_l\)的乘积了。同理，第16行，\(\mathbf{c^T-c_B^TB^{-1}A}\)，因为也是只有\(c_e\)不为0，所以就和他对应的\(\mathbf{A}\)的第\(l\)行相乘了。

到此，终于介绍完了单纯形算法。其余还有一些要注意的地方，好比必定要注意检验数和原目标函数的\(\mathbf{c}\)是彻底不同的概念，在原约束为不等式，须要加松弛变量的状况下，他们可能相等，但内心必定要区分它们，同时，这种状况下，基很容易找，就是松弛变量的那几列构成的单位阵。可是若是原约束是等式，就须要本身找基，而且这时检验数每每就和目标函数参数不一样了。

最后，本文所用的截图均来自中科院计算所B老师的课程PPT，本人在学习该课程时也受益良多，对单纯形算法也钻研了比较长的时间，所以撰写本文，但愿给你们一个学习参考！其中可能有错误之处，欢迎指正讨论。

对于原创博文：如需转载请注明出处http://www.cnblogs.com/Kenneth-Wong/