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

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

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


写作背景

大家都知道ML的基础是统计和优化,作者本人这学期选了最优化计算(还没开课),顺道旁听了一门组合优化,这次就写写组合优化的入门级算法---单纯形算法,虽然这不是一个多项式算法,但是它的应用面极为广泛。本来没想投入这么多精力进去的,然后自己就沦陷了,都快抛弃主课了。数学专业上课讲的东西更偏算法设计的思想和主要流程,本人算是半个数学爱好者,除了学习数学的思维方式,也学一些CS知识敲代码,所以追求一点细节以及代码实现,写出来分享一些我遇到的问题。

虽然本人不是一个数学专业的学生,但是这篇文章还是需要一些数学知识储备的,同时也不准备举例子,打算讲清楚算法原理和推导,害怕推导和公式者慎入。

1949年:G.B.Dantzig提出了单纯形方法。关于历史考究什么的不想多提,可以自行搜索。


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

线性规划 (LP) 问题的标准型:

min\ z = c^T x\\ s.t.\begin{cases} Ax=b\\ x\ge0 \end{cases}

其中 x\in R^n, c\in R^n, b\in R^m, A\in R^{m \times n}

一般总假设 b\ge0,\ (A,b,c) 的元素都为整数(线性规划多应用在实际问题,所以要求为整数也可以理解), rank(A)=m

可行域: P = \{ x\in R^n : Ax=b, x\ge0\}

最优解集: P^*=\{x\in P:x(LP) 的最优解 \}


二 定理和推论

如果对仿射集,凸集,凸包和多胞形这些概念有一定的了解的话,看一些证明也有助于理解。如果对定理证明不理解也没什么关系,做为一名工科生,能够理解大致意思和应用就是最好的结果。

定理1x 是可行域 P 的一个顶点 \Leftrightarrow x 的正分量对应的 A 中的各列是线性独立的。

  • 对矩阵 A 进行分解(可比对高斯消去法解方程组时,寻找主元变量的过程),移列使得前 m列线性不相关A=\begin{bmatrix} B, N\end{bmatrix} , B\in R^{m \times m} 为满秩矩阵,理解为矩阵 A 空间的一组基向量x=\begin{bmatrix} x_B\\ x_N \end{bmatrix} 。则有:

\begin{equation} Ax=b\ \Leftrightarrow \ Bx_B+Nx_N=b \ \Leftrightarrow \ x_B+B^{-1}Nx_N=B^{-1}b \end{equation}

\Leftrightarrow \ x_B=B^{-1}b-B^{-1}Nx_N ,令 x_N=0 ,若 x= \begin{bmatrix} B^{-1}b\\0 \end{bmatrix} \ge 0 ,称之为 (LP) 的一个基可行解。

  • B 满秩, A=\begin{bmatrix} B, N\end{bmatrix} 中有 m 个基向量组成线性无关组, x=\begin{bmatrix} x_B\\ x_N \end{bmatrix} 中有 m 个主元变量, n-m 个自由变量。通过上述分解可有以下推论:

推论1x(LP) 的一个基可行解 \Leftrightarrow x \in PP 的一个顶点。

推论2(LP) 的可行域至多有 C^m_n 个顶点。

如果 x_B >0,\ x_N=0 ,则称一个基可行解 x=\begin{bmatrix} x_B\\ x_N \end{bmatrix} 非退化。

由于可能有退化的情况,所以基可行解和顶点不是一一映射:任给一个基可行解,存在唯一的一个顶点与之对应;对于 P 的一个顶点,可能有多个基可行解与之对应。

如果它所对应的所有基可行解都是非退化的,则称这个线性问题非退化。

考虑只有一个基变量不同的基可行解,可知这两个基可行解对应的顶点相邻。


三 单纯形算法框架

给定一个非退化的基可行解 \bar x=\begin{bmatrix} x_B\\ x_N \end{bmatrix}注意这里没有给出求基可行解方法,后面会说明怎么求初始迭代的基可行解),对应的可行基为 B 则等式约束 Ax=b 变为:

x_B+B^{-1}Nx_N=B^{-1}b\\ x_B=B^{-1}b-B^{-1}Nx_N

目标函数

\begin{align} c^Tx &= c_B^Tx_B+c_N^Tx_N\\ &=c_B^T(B^{-1}b-B^{-1}Nx_N)+c_N^Tx_N\\ &= c_B^TB^{-1}b-(c_B^TB^{-1}N-c_N^T)x_N \\ \end{align}

规划等价于

min\ c_B^TB^{-1}b-(c_B^TB^{-1}N-c_N^T)x_N\\ s.t.\begin{cases} x_B=B^{-1}b-B^{-1}Nx_N\\ x\ge0 \end{cases}

做以下记号:

B^{-1}b= \begin{bmatrix} \alpha_1\\ \vdots\\ \alpha_m \end{bmatrix}B^{-1}N= \begin{bmatrix} \beta_{1,m+1} &\beta_{1,m+2} &\cdots &\beta_{1,n} \\ \beta_{2,m+1} &\beta_{2,m+2} &\cdots &\beta_{2,n}\\ \vdots & \vdots & \ddots & \vdots \\ \beta_{m,m+1} &\beta_{m,m+2} &\cdots &\beta_{m,n} \end{bmatrix}

\begin{bmatrix} \lambda_{m+1} &\lambda_{m+2} &\cdots &\lambda_n \end{bmatrix} =c_B^TB^{-1}N-c_N^T ,\ f_0=c_B^TB^{-1}b

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

min\ f_0-(\lambda_{m+1}x_{m+1}+\lambda_{m+2}x_{m+2}+\cdots+\lambda_nx_n)\\ s.t.\begin{cases} x_1=\alpha_1-(\beta_{1,m+1}x_{m+1}+\beta_{1,m+2}x_{m+2}+\cdots +\beta_{1,n}x_n\\ x_2=\alpha_2-(\beta_{2,m+1}x_{m+1}+\beta_{2,m+2}x_{m+2}+\cdots +\beta_{2,n}x_n\\ \cdots \cdots\\ x_m=\alpha_m-(\beta_{m,m+1}x_{m+1}+\beta_{m,m+2}x_{m+2}+\cdots +\beta_{m,n}x_n\\ x_i\ge0 (i=1,2,\cdots,n) \end{cases}

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

x_N= \begin{bmatrix} x_{m+1}\\ \vdots\\ x_n \end{bmatrix}=0,对应的基可行解为 \begin{bmatrix} B^{-1}b\\ 0 \end{bmatrix} ,目标函数值为 f_0

  1. 在当前给定基可行解 \bar x 的情况下,若有一个 \lambda_i >0 ,不妨设 \lambda_{m+1} > 0 ,则可令 x_{m+1} 从0上升到某个 \theta>0 ,使得目标函数值下降。说明当前解非最优解。
  2. 如果检验数 \lambda_{m+1} \le 0,\lambda_{m+2}\le0,\cdots,\lambda_{m+n}\le0,则基 \begin{bmatrix} x_1,x_2,\cdots,x_m \end{bmatrix} 对应的基可行解 x_0 = \begin{bmatrix} \alpha_1,\alpha_2,\cdots,\alpha_m,0,\cdots,0 \end{bmatrix}^T 是最优解。
  3. 如果目标函数有下界且存在一个检验数 \lambda_{m+k} >0 ,则非基变量 x_{m+k} 对应的系数 \beta_{1,m+k},\beta_{2,m+k},\cdots,\beta_{m,m+k} 中至少有一个大于0.

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


四 细节解释

  • 求基可行解的方法

最开始听完课就在想这怎么划分 A=\begin{bmatrix} B, N\end{bmatrix} ,如何找一组线性无关的向量,进行高斯消去找线性无关组?可是程序化怎么实现呢?

后来去搜了资料,简单易行的方法就是在标准型上构造人工向量。

原始问题标准型
对约束构造人工变量

聪明的你们可能会发现约束已经变了,那目标问题求解也就变了。没错,所以两种通用的方法就出来了---大M法和两阶段法。我只简单介绍一种大M法,很容易看懂,两阶段法可以去搜资料看教材。

M是一个充分大的数,目标函数求最大就取减号,求最小就取加号。在这个问题中要有最优解,人工变量就必须取0,否则,原问题无解。这样就很顺利的得到了一个基可行解,然后代入单纯形算法进行迭代

  • 入基变量的选取:

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

  • 出基变量的选取:

最小比值法选取出基变量,当选取完入基变量 \lambda_{m+k} 后,取 x_{m+k}=min \ \frac{\alpha_i}{\beta_{i,m+k}} ,相应的 x_i 做为出基变量置为0,取最小的原则是可以保证 x_i\ge0 (i=1,2,\cdots,n) ,防止出现非负变量。


五 结语

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



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



欢迎原链接转发,转载请前往@zhihu.com/people/leilei的主页获取信息,盗版必究。


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


如果你是运筹学/人工智能硕博或在读,请在下图的公众号后台留言:“加微信群”。系统会自动辨认你的关键字,并提示您进一步的加群要求和步骤,邀请您进全球运筹或AI学者群(群内学界、业界大佬云集)。

同时我们有:【运筹学|优化爱好者】【供应链|物流】【人工智能】【数据科学|分析】千人QQ群,想入群的小伙伴可以关注下方公众号点击“加入社区”按钮,获得入群传送门。

学术界|工业界招聘、征稿等信息免费发布,请见下图:

编辑于 2018-06-08

文章被以下专栏收录