论坛 › 技术讨论 Technology Answers › 【学界】单纯形算法的原理和示例实现

ainatipen 发表于 2021-12-16 20:36

【学界】单纯形算法的原理和示例实现

作者简介： @磊磊 本科河北工业大学电子科学与技术专业，现在中科院半导体研究所读硕士，目前研一，方向NIRS或CV相关。自己本着兴趣学习CS,ML和统计优化相关知识，写一些自己的学习笔记为正在学习路上的小伙伴做做铺垫，内容可能偏数学一点，因为本人算是半个数学爱好者，打的math基础也更多是用在算法理解方面。
欢迎原链接转发，付费转载请前往 @留德华叫兽 的主页获取信息，盗版必究。
敬请关注和扩散本专栏及同名公众号，会邀请全球知名学者陆续发布运筹学、人工智能中优化理论等相关干货、知乎Live及行业动态：
『运筹OR帷幄』大数据人工智能时代的运筹学--知乎专栏
写作背景

大家都知道ML的基础是统计和优化，作者本人这学期选了最优化计算(还没开课)，顺道旁听了一门组合优化，这次就写写组合优化的入门级算法---单纯形算法，虽然这不是一个多项式算法，但是它的应用面极为广泛。本来没想投入这么多精力进去的，然后自己就沦陷了，都快抛弃主课了。数学专业上课讲的东西更偏算法设计的思想和主要流程，本人算是半个数学爱好者，除了学习数学的思维方式，也学一些CS知识敲代码，所以追求一点细节以及代码实现，写出来分享一些我遇到的问题。
虽然本人不是一个数学专业的学生，但是这篇文章还是需要一些数学知识储备的，同时也不准备举例子，打算讲清楚算法原理和推导，害怕推导和公式者慎入。
1949年：G.B.Dantzig提出了单纯形方法。关于历史考究什么的不想多提，可以自行搜索。

一线性规划问题的标准形式

线性规划问题的标准型：

https://www.zhihu.com/equation?tex=min%5C+z+%3D+c%5ET+x%5C%5C+s.t.%5Cbegin%7Bcases%7D+Ax%3Db%5C%5C+x%5Cge0+%5Cend%7Bcases%7D
其中 https://www.zhihu.com/equation?tex=x%5Cin+R%5En%2C+c%5Cin+R%5En%2C+b%5Cin+R%5Em%2C+A%5Cin+R%5E%7Bm+%5Ctimes+n%7D
一般总假设 https://www.zhihu.com/equation?tex=b%5Cge0%2C%5C+%28A%2Cb%2Cc%29 的元素都为整数（线性规划多应用在实际问题，所以要求为整数也可以理解）， https://www.zhihu.com/equation?tex=rank%28A%29%3Dm
可行域： https://www.zhihu.com/equation?tex=P+%3D+%5C%7B+x%5Cin+R%5En+%3A+Ax%3Db%2C+x%5Cge0%5C%7D
最优解集： https://www.zhihu.com/equation?tex=P%5E%2A%3D%5C%7Bx%5Cin+P%3Ax 是的最优解 https://www.zhihu.com/equation?tex=%5C%7D

二定理和推论

如果对仿射集，凸集，凸包和多胞形这些概念有一定的了解的话，看一些证明也有助于理解。如果对定理证明不理解也没什么关系，做为一名工科生，能够理解大致意思和应用就是最好的结果。
定理1：是可行域的一个顶点 https://www.zhihu.com/equation?tex=%5CLeftrightarrow+x 的正分量对应的中的各列是线性独立的。





[*]对矩阵进行分解（可比对高斯消去法解方程组时，寻找主元变量的过程），移列使得前 列线性不相关，https://www.zhihu.com/equation?tex=A%3D%5Cbegin%7Bbmatrix%7D+B%2C+N%5Cend%7Bbmatrix%7D+%2C+B%5Cin+R%5E%7Bm+%5Ctimes+m%7D 为满秩矩阵，理解为矩阵空间的一组基向量，。则有：

https://www.zhihu.com/equation?tex=%5Cbegin%7Bequation%7D+Ax%3Db%5C+%5CLeftrightarrow+%5C+Bx_B%2BNx_N%3Db+%5C+%5CLeftrightarrow+%5C+x_B%2BB%5E%7B-1%7DNx_N%3DB%5E%7B-1%7Db+%5Cend%7Bequation%7D

https://www.zhihu.com/equation?tex=%5CLeftrightarrow+%5C+x_B%3DB%5E%7B-1%7Db-B%5E%7B-1%7DNx_N ，令 https://www.zhihu.com/equation?tex=x_N%3D0 ，若 https://www.zhihu.com/equation?tex=x%3D+%5Cbegin%7Bbmatrix%7D+B%5E%7B-1%7Db%5C%5C0+%5Cend%7Bbmatrix%7D+%5Cge+0+ ，称之为的一个基可行解。

[*] 满秩，中有个基向量组成线性无关组，中有个主元变量， https://www.zhihu.com/equation?tex=n-m 个自由变量。通过上述分解可有以下推论：
推论1：是 https://www.zhihu.com/equation?tex=%EF%BC%88LP%29 的一个基可行解 https://www.zhihu.com/equation?tex=%5CLeftrightarrow+x+%5Cin+P 是的一个顶点。
推论2：的可行域至多有 https://www.zhihu.com/equation?tex=C%5Em_n 个顶点。
如果 https://www.zhihu.com/equation?tex=x_B+%3E0%2C%5C+x_N%3D0 ，则称一个基可行解非退化。
由于可能有退化的情况，所以基可行解和顶点不是一一映射：任给一个基可行解，存在唯一的一个顶点与之对应；对于的一个顶点，可能有多个基可行解与之对应。
如果它所对应的所有基可行解都是非退化的，则称这个线性问题非退化。
考虑只有一个基变量不同的基可行解，可知这两个基可行解对应的顶点相邻。

三单纯形算法框架

给定一个非退化的基可行解 https://www.zhihu.com/equation?tex=%5Cbar+x%3D%5Cbegin%7Bbmatrix%7D+x_B%5C%5C+x_N+%5Cend%7Bbmatrix%7D （注意这里没有给出求基可行解方法，后面会说明怎么求初始迭代的基可行解）,对应的可行基为则等式约束 https://www.zhihu.com/equation?tex=Ax%3Db 变为：

https://www.zhihu.com/equation?tex=x_B%2BB%5E%7B-1%7DNx_N%3DB%5E%7B-1%7Db%5C%5C+x_B%3DB%5E%7B-1%7Db-B%5E%7B-1%7DNx_N
目标函数

https://www.zhihu.com/equation?tex=%5Cbegin%7Balign%7D+c%5ETx+%26%3D+c_B%5ETx_B%2Bc_N%5ETx_N%5C%5C+%26%3Dc_B%5ET%28B%5E%7B-1%7Db-B%5E%7B-1%7DNx_N%29%2Bc_N%5ETx_N%5C%5C+%26%3D+c_B%5ETB%5E%7B-1%7Db-%28c_B%5ETB%5E%7B-1%7DN-c_N%5ET%29x_N+%5C%5C+%5Cend%7Balign%7D
规划等价于

https://www.zhihu.com/equation?tex=min%5C+c_B%5ETB%5E%7B-1%7Db-%28c_B%5ETB%5E%7B-1%7DN-c_N%5ET%29x_N%5C%5C+s.t.%5Cbegin%7Bcases%7D+x_B%3DB%5E%7B-1%7Db-B%5E%7B-1%7DNx_N%5C%5C+x%5Cge0+%5Cend%7Bcases%7D
做以下记号：

https://www.zhihu.com/equation?tex=B%5E%7B-1%7Db%3D+%5Cbegin%7Bbmatrix%7D+%5Calpha_1%5C%5C+%5Cvdots%5C%5C+%5Calpha_m+%5Cend%7Bbmatrix%7D ， https://www.zhihu.com/equation?tex=B%5E%7B-1%7DN%3D+%5Cbegin%7Bbmatrix%7D+%5Cbeta_%7B1%2Cm%2B1%7D+%26%5Cbeta_%7B1%2Cm%2B2%7D+%26%5Ccdots+%26%5Cbeta_%7B1%2Cn%7D+%5C%5C+%5Cbeta_%7B2%2Cm%2B1%7D+%26%5Cbeta_%7B2%2Cm%2B2%7D+%26%5Ccdots+%26%5Cbeta_%7B2%2Cn%7D%5C%5C+%5Cvdots+%26+%5Cvdots+%26+%5Cddots+%26+%5Cvdots+%5C%5C+%5Cbeta_%7Bm%2Cm%2B1%7D+%26%5Cbeta_%7Bm%2Cm%2B2%7D+%26%5Ccdots+%26%5Cbeta_%7Bm%2Cn%7D+%5Cend%7Bbmatrix%7D ，

https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+%5Clambda_%7Bm%2B1%7D+%26%5Clambda_%7Bm%2B2%7D+%26%5Ccdots+%26%5Clambda_n+%5Cend%7Bbmatrix%7D+%3Dc_B%5ETB%5E%7B-1%7DN-c_N%5ET+%2C%5C+f_0%3Dc_B%5ETB%5E%7B-1%7Db
则上述线性规划问题就变成：

https://www.zhihu.com/equation?tex=min%5C+f_0-%28%5Clambda_%7Bm%2B1%7Dx_%7Bm%2B1%7D%2B%5Clambda_%7Bm%2B2%7Dx_%7Bm%2B2%7D%2B%5Ccdots%2B%5Clambda_nx_n%29%5C%5C+s.t.%5Cbegin%7Bcases%7D+x_1%3D%5Calpha_1-%28%5Cbeta_%7B1%2Cm%2B1%7Dx_%7Bm%2B1%7D%2B%5Cbeta_%7B1%2Cm%2B2%7Dx_%7Bm%2B2%7D%2B%5Ccdots+%2B%5Cbeta_%7B1%2Cn%7Dx_n%5C%5C+x_2%3D%5Calpha_2-%28%5Cbeta_%7B2%2Cm%2B1%7Dx_%7Bm%2B1%7D%2B%5Cbeta_%7B2%2Cm%2B2%7Dx_%7Bm%2B2%7D%2B%5Ccdots+%2B%5Cbeta_%7B2%2Cn%7Dx_n%5C%5C+%5Ccdots+%5Ccdots%5C%5C+x_m%3D%5Calpha_m-%28%5Cbeta_%7Bm%2Cm%2B1%7Dx_%7Bm%2B1%7D%2B%5Cbeta_%7Bm%2Cm%2B2%7Dx_%7Bm%2B2%7D%2B%5Ccdots+%2B%5Cbeta_%7Bm%2Cn%7Dx_n%5C%5C+x_i%5Cge0+%28i%3D1%2C2%2C%5Ccdots%2Cn%29+%5Cend%7Bcases%7D
学数学的人为了表示方便，一般把上述问题转换为一个单纯形表去表示（其实就是写成矩阵形式），但是我觉得就这样看容易分析算法，然后去用代码迭代实现了。毕竟本人不是专门学数学的，怎么方便理解怎么来，图也顺道贴上，有点丑别见怪。


当 https://www.zhihu.com/equation?tex=x_N%3D+%5Cbegin%7Bbmatrix%7D+x_%7Bm%2B1%7D%5C%5C+%5Cvdots%5C%5C+x_n+%5Cend%7Bbmatrix%7D%3D0，对应的基可行解为 https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+B%5E%7B-1%7Db%5C%5C+0+%5Cend%7Bbmatrix%7D ，目标函数值为 https://www.zhihu.com/equation?tex=f_0 。

[*]在当前给定基可行解 https://www.zhihu.com/equation?tex=%5Cbar+x 的情况下，若有一个 https://www.zhihu.com/equation?tex=%5Clambda_i+%3E0 ，不妨设 https://www.zhihu.com/equation?tex=%5Clambda_%7Bm%2B1%7D+%3E+0 ，则可令 https://www.zhihu.com/equation?tex=x_%7Bm%2B1%7D 从0上升到某个 https://www.zhihu.com/equation?tex=%5Ctheta%3E0 ，使得目标函数值下降。说明当前解非最优解。
[*]如果检验数 https://www.zhihu.com/equation?tex=%5Clambda_%7Bm%2B1%7D+%5Cle+0%2C%5Clambda_%7Bm%2B2%7D%5Cle0%2C%5Ccdots%2C%5Clambda_%7Bm%2Bn%7D%5Cle0,则基 https://www.zhihu.com/equation?tex=%5Cbegin%7Bbmatrix%7D+x_1%2Cx_2%2C%5Ccdots%2Cx_m+%5Cend%7Bbmatrix%7D 对应的基可行解 https://www.zhihu.com/equation?tex=x_0+%3D+%5Cbegin%7Bbmatrix%7D+%5Calpha_1%2C%5Calpha_2%2C%5Ccdots%2C%5Calpha_m%2C0%2C%5Ccdots%2C0+%5Cend%7Bbmatrix%7D%5ET 是最优解。
[*]如果目标函数有下界且存在一个检验数 https://www.zhihu.com/equation?tex=%5Clambda_%7Bm%2Bk%7D+%3E0 ，则非基变量 https://www.zhihu.com/equation?tex=x_%7Bm%2Bk%7D 对应的系数 https://www.zhihu.com/equation?tex=%5Cbeta_%7B1%2Cm%2Bk%7D%2C%5Cbeta_%7B2%2Cm%2Bk%7D%2C%5Ccdots%2C%5Cbeta_%7Bm%2Cm%2Bk%7D 中至少有一个大于0.
（通俗一点就是如果当前解不是最优解，寻找入基变量和出基变量，即寻找其相邻的顶点，在这里就需要注意一点，为什么要求问题非退化，退化的情况有多个基可行解对应一个顶点，迭代更新时会陷入循环出不去）

四细节解释


[*]求基可行解的方法：
最开始听完课就在想这怎么划分，如何找一组线性无关的向量，进行高斯消去找线性无关组？可是程序化怎么实现呢？
后来去搜了资料，简单易行的方法就是在标准型上构造人工向量。



原始问题标准型



对约束构造人工变量

聪明的你们可能会发现约束已经变了，那目标问题求解也就变了。没错，所以两种通用的方法就出来了---大M法和两阶段法。我只简单介绍一种大M法，很容易看懂，两阶段法可以去搜资料看教材。
M是一个充分大的数，目标函数求最大就取减号，求最小就取加号。在这个问题中要有最优解，人工变量就必须取0，否则，原问题无解。这样就很顺利的得到了一个基可行解，然后代入单纯形算法进行迭代。

[*]入基变量的选取：
最简单的选取规则就是：目标函数求max,检验数大的为入基变量,目标函数求min,检验数小的为入基变量。当然还有别的选取方式，可以去搜搜相关的paper去看看。这里只是一个简单的科普。

[*]出基变量的选取：
最小比值法选取出基变量，当选取完入基变量 https://www.zhihu.com/equation?tex=%5Clambda_%7Bm%2Bk%7D 后，取 https://www.zhihu.com/equation?tex=x_%7Bm%2Bk%7D%3Dmin+%5C+%5Cfrac%7B%5Calpha_i%7D%7B%5Cbeta_%7Bi%2Cm%2Bk%7D%7D ，相应的 https://www.zhihu.com/equation?tex=x_i 做为出基变量置为0，取最小的原则是可以保证 https://www.zhihu.com/equation?tex=x_i%5Cge0+%28i%3D1%2C2%2C%5Ccdots%2Cn%29 ，防止出现非负变量。

五 结语

这篇文章只是写了单纯形算法的设计框架和原理，并没有完全概括出来设计思想，算是入门型教学，后面我还会更新一篇文章，结合对偶单纯形和原始对偶算法来提炼出整体的设计思想，希望到时候读者可以看得更清楚些，这篇是后面的基础和铺垫。

MATLAB code Demo
%匆忙之间写的代码，我也没有测试有没有问题，就是提供一个示例
%大M法构造人工变量的个数可以考虑根据A有没有单位列进行减少这个可以添加上

function = BigMsimplexmin(A, c, b)
%A   m*n
%x   n*1
%b   m*1
%c   n*1

%大M法构造基可行解
M=5000;
= size(A);
A = ;
c = ;
x = ;

while (1)
    %分离基变量和非基变量
    index_B = find(x ~=0);
    index_N = find(x == 0);
    B = A(:, index_B);
    N = A(:, index_N);
%   x_B = x(index_B);
    x_N = x(index_N);
    c_B = c(index_B);
    c_N = c(index_N);
   
    %计算lambda，alpha，beta，f
    lambda = (c_B'/B)*N-c_N';    %1*(n-m)
    alpha = B\b;
    beta = B\N;
    f = c_B'*alpha;
   
    %lambda都小于0, 当前解为最优解,返回x和f
    if isempty(find(sign(lambda) > 0, 1))
      x = x;
      return
    end
   
    % 选最大的正检验数作为进基变量
    [~, k] = max(lambda);
   
    %更新入基变量
    x_N(k) = min(alpha./beta(:, k));
   
    %更新基变量,同时置出基变量为0
    x_B = alpha - beta(:, k).*x_N;
   
    %新的x0
    x = ;
end

end

欢迎原链接转发，转载请前往@https://www.zhihu.com/people/leilei-1996/posts的主页获取信息，盗版必究。

以上『运筹OR帷幄』专栏所有文章都会同步发送至留德华叫兽的头条主页， 以及同名微信公众号，目前预计受众10w +

如果你是运筹学/人工智能硕博或在读，请在下图的公众号后台留言：“加微信群”。系统会自动辨认你的关键字，并提示您进一步的加群要求和步骤，邀请您进全球运筹或AI学者群（群内学界、业界大佬云集）。
同时我们有：【运筹学|优化爱好者】【供应链|物流】【人工智能】【数据科学|分析】千人QQ群，想入群的小伙伴可以关注下方公众号点击“加入社区”按钮，获得入群传送门。
学术界|工业界招聘、征稿等信息免费发布，请见下图：

JamesB 发表于 2021-12-16 20:36

写的挺好的 附代码很赞

unityloverz 发表于 2021-12-16 20:38

代码不敢保证，今天匆忙写出来了，代码考虑不是很完善，我还没测试过，就是提供个demo。

kyuskoj 发表于 2021-12-16 20:44

demo已经足够啦

IT圈老男孩1 发表于 2021-12-16 20:49

感谢投稿～

RhinoFreak 发表于 2021-12-16 20:56

半导体所也开始搞CV了…

stonstad 发表于 2021-12-16 20:59

很少，我们组就是一直有这个方向，还是当时王守觉先生在的时候就有了

zt3ff3n 发表于 2021-12-16 21:05

求问 ，如果B不可逆 如何处理？

fwalker 发表于 2021-12-16 21:10

找的就是可逆的B

Arzie100 发表于 2021-12-16 21:10

您好，想问一下为什么一个检测数lamda大于0那么对应的非基变量beta系数大于0呢

查看完整版本: 【学界】单纯形算法的原理和示例实现