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

ainatipen

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

写作背景

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

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

线性规划  问题的标准型：


其中
一般总假设 的元素都为整数（线性规划多应用在实际问题，所以要求为整数也可以理解），
可行域：
最优解集： 是  的最优解

二  定理和推论

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

• 对矩阵  进行分解（可比对高斯消去法解方程组时，寻找主元变量的过程），移列使得前 列线性不相关， 为满秩矩阵，理解为矩阵  空间的一组基向量，  。则有：
  

，令 ，若 ，称之为  的一个基可行解。

• 满秩，  中有  个基向量组成线性无关组，  中有  个主元变量， 个自由变量。通过上述分解可有以下推论：
  

推论1：  是 的一个基可行解 是  的一个顶点。
推论2：  的可行域至多有 个顶点。
如果 ，则称一个基可行解  非退化。
由于可能有退化的情况，所以基可行解和顶点不是一一映射：任给一个基可行解，存在唯一的一个顶点与之对应；对于  的一个顶点，可能有多个基可行解与之对应。
如果它所对应的所有基可行解都是非退化的，则称这个线性问题非退化。
考虑只有一个基变量不同的基可行解，可知这两个基可行解对应的顶点相邻。

三  单纯形算法框架

给定一个非退化的基可行解 （注意这里没有给出求基可行解方法，后面会说明怎么求初始迭代的基可行解）,对应的可行基为  则等式约束 变为：


目标函数


规划等价于


做以下记号：

， ，


则上述线性规划问题就变成：


学数学的人为了表示方便，一般把上述问题转换为一个单纯形表去表示（其实就是写成矩阵形式），但是我觉得就这样看容易分析算法，然后去用代码迭代实现了。毕竟本人不是专门学数学的，怎么方便理解怎么来，图也顺道贴上，有点丑别见怪。


当 ，对应的基可行解为 ，目标函数值为 。

• 在当前给定基可行解 的情况下，若有一个 ，不妨设 ，则可令 从0上升到某个 ，使得目标函数值下降。说明当前解非最优解。
  
• 如果检验数 ,则基 对应的基可行解 是最优解。
  
• 如果目标函数有下界且存在一个检验数 ，则非基变量 对应的系数 中至少有一个大于0.
  

（通俗一点就是如果当前解不是最优解，寻找入基变量和出基变量，即寻找其相邻的顶点，在这里就需要注意一点，为什么要求问题非退化，退化的情况有多个基可行解对应一个顶点，迭代更新时会陷入循环出不去）

四  细节解释

• 求基可行解的方法：
  

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



原始问题标准型



对约束构造人工变量

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

• 入基变量的选取：
  

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

• 出基变量的选取：
  

最小比值法选取出基变量，当选取完入基变量 后，取 ，相应的 做为出基变量置为0，取最小的原则是可以保证 ，防止出现非负变量。

五 结语

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

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

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

%大M法构造基可行解
M=5000;
[m,n] = size(A);
A = [A, eye(m)];
c = [c;repmat(M,m,1)];
x = [zeros(n,1);b];

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 = [x_B; x_N];
end

end

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

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

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

JamesB

写的挺好的 附代码很赞

unityloverz

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

kyuskoj

demo已经足够啦

IT圈老男孩1

感谢投稿～

RhinoFreak

半导体所也开始搞CV了…

stonstad

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

zt3ff3n

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

fwalker

找的就是可逆的B

Arzie100

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