【学界】原始-对偶算法的设计原理

【学界】原始-对偶算法的设计原理

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


这次就接着上一篇磊磊:单纯形算法的原理和示例实现继续写,建议读这篇文章前,最好看看前面的单纯形那篇文章。

这两篇文章的信息量有点大,我希望有心学习的小伙伴可以自己推导推导加深理解。


1. 拉格朗日对偶

如果大家对本科学的高等数学还有印象的话,应该记得讲过拉格朗日乘子法求极值,一般流程是把约束问题写成拉格朗日乘子形式,然后求导求极值,这个东西很好用,我自己在学习机器学习过程中有时也用拉格朗日乘子法进行推导算法,学SVM避不开这个东西,在推导PCA和

磊磊:高斯混合模型(GMM)的两种详解及简化时同样可以用到。这里提出来是要把原始问题转换为对偶问题。

LP问题的标准型:

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

采用拉格朗日乘子法 \lambda \in R^m ,对偶函数为

\begin{align}d(\lambda)&=\min \limits_{x\ge0}\{ c^Tx+\lambda^T(b-Ax) \} =b^T\lambda+\min\limits_{x\ge0}(c-A^T\lambda)^Tx\\ &=\begin{cases}b^T\lambda,\ \ \mbox{if}\ A^T\lambda\leq c\\ -\infty, \ \ \mbox{othesize} \end{cases} \end{align}

所以就引出了对偶问题(LD)

\max \ b^T\lambda \\ s.t.\begin{cases} A^T\lambda \leq c\\ \lambda \ \ \ \mbox{no constraint } \end{cases}

有了LP标准型和它的LD形式,对于非标准型的形式可以引入松弛变量变换为标准型,然后套用上面的形式,这一点在一些推导或者变换问题时非常有用,下面举个例子展示一下如何引入松弛变量:

假设原始问题为:

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

在目标求值和约束方程中引入松弛变量 s \ge 0 可得 z=\begin{bmatrix} c^T, 0 \end{bmatrix} \begin{bmatrix} x\\ s \end{bmatrix} \begin{bmatrix} A, -I \end{bmatrix} \begin{bmatrix} x\\ s \end{bmatrix} =b ,转换为对偶问题时有约束 \begin{bmatrix} A, -I \end{bmatrix}^T\lambda \leq \begin{bmatrix} c\\0 \end{bmatrix} ,则对偶问题为:

\max \ b^T\lambda \\ s.t.\begin{cases} A^T\lambda \leq c\\ \lambda \ge0 \end{cases}

遇到带有不等式的情况就算拉格朗日乘子法如何转换没记住,只记住一个标准型的原始对偶转换,再利用引入松弛变量也可以做到游刃有余。


2. 对偶性质

弱对偶性:

如果 \bar x_i(i=1,2,\cdots,n) 是原问题可行解, \bar y_j(j=1,2,\cdots,m) 是其对偶问题可行解,则恒有 c^Tx \ge b^Ty .

(LP)\min \ c^Tx \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (LD)\max \ b^Ty\\ s.t.\begin{cases} Ax=b\\ x\ge0 \end{cases} \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ s.t.\begin{cases} A^Ty \leq c\\ y \ \ \ \mbox{no constraint } \end{cases}

证明:

\sum_{i=1}^{n}{c_i\bar x_i} \ge \sum_{i=1}^{n}{\sum_{j=1}^{m}{a_{ij} \bar y_j \bar x_i}} \ , \ \sum_{j=1}^{m}{b_j\bar y_j} = \sum_{j=1}^{m}{\sum_{i=1}^{n}{a_{ij} \bar x_i \bar y_j}} \ 可证。

由弱对偶性,可得出以下推论

  • 原问题任意可行解的目标函数值是其对偶问题目标函数值的上界,反之,对偶问题任意可行解的目标函数值是其原问题目标函数值的下界。
  • 在一对对偶问题(LP)和(LD)中,若其中一个问题可行但目标函数无界,则另一个问题无可行解。
  • 在一对对偶问题(LP)和(LD)中,若一个(如LP)可行,而另一个(如LD)不可行,则该可行的问题目标函数值无界。


最优性:

如果 x 是原问题的可行解, y 是对偶问题的可行解,且 c^Tx=b^Ty ,则 x,y 分别是原始问题和对偶问题的最优解。

证:设 x^*,y^*分别是原始问题和对偶问题的最优解,则有 c^Tx\ge c^Tx^*\ge b^Ty^*\ge b^Ty

又有 c^Tx=b^Ty ,则 c^Tx=c^Tx^*=b^Ty^*=b^Tyx,y 分别是原始问题和对偶问题的最优解。


强对偶性:

若原问题及其对偶问题均具有可行解,则两者均具有最优解,且他们最优解的目标函数值相等。

证:由于两者均有可行解,根据弱对偶性的推论,对原问题的目标函数值具有上界,对偶问题的目标函数值具有下界,因此两者均具有最优解。参考第4节的再看单纯形算法,假设原问题的最优解为 x ,则必然满足检验数 \lambda = c_B^TB^{-1}N-c_N^T \le0 ,记 y^T=c_B^TB^{-1} ,有 A^Ty-c=\begin{bmatrix} 0\\ c_B^TB^{-1}N-c_N^T \end{bmatrix}\le 0 , A^Ty \le cy 是其对偶问题的可行解,此时对偶问题的 b^Ty = y^Tb = c_B^TB^{-1}b=c^Tx ,由最优性知,这时两者的解均为最优解。


3. 互补松弛性条件

假设 x,y 分别为原始和对偶问题的最优解在对偶问题中引入松弛变量 s ,可得如下方程组:

\begin{cases} Ax-b=0 \ \ \mbox{LP可行域}\\ A^Ty+s-c=0 \ \ \mbox{LD可行域}\\ c^Tx-b^Ty=0 \ \ \mbox{强对偶性}\\ x \ge 0, \ s\ge 0 \end{cases}

c^Tx=x^T(A^Ty+s)=x^TA^Ty+x^Ts=b^Ty+x^Ts \Leftrightarrow x^Ts=0\\ \Leftrightarrow \sum_{i=1}^{n}{x_is_i} \Leftrightarrow x_is_i=0(i=1,2 ,\cdots,n)\Leftrightarrow x_i(c_i-A_i^Ty)=0

x_i(c_i-A_i^Ty)=0 即是互补松弛性条件。这个东西还有一个用处,当我们求出一个问题的最优解和最优目标函数值时,可以通过对偶松弛性条件和强对偶性去求其对偶问题的最优解,在后面会提到有应用。

这一节的内容别看少,但是这是算法设计的核心思想。


4. 再看单纯形算法

看过前一篇文章的话,应该会记得我们在进行单纯形算法设计时的做法。记不住的话请回看,这里可能会简写一些。

  • 在基可行解的情况下要求x=\begin{bmatrix} x_B,\\x_N \end{bmatrix}=\begin{bmatrix} B^{-1}b\\0\end{bmatrix} \ge0说明原始问题可行
  • 目标函数值 \min \ f_0=c_B^TB^{-1}b ,记 y^T=c_B^TB^{-1} , 则有c^Tx=c_B^TB^{-1}b+0=b^Ty说明满足对偶间隙为0。
  • y^TA-c^T=\begin{bmatrix} 0,\ c_B^TB^{-1}N-c_N^T \end{bmatrix}是不是发现了 c_B^TB^{-1}N-c_N^T 就是单纯形算法中要改善的检验数 \lambda ,改善对偶不可行。

对比上一节互补性松弛条件下面的方程组,单纯形算法的设计就是满足了三个条件中的两个,去改善另一个条件。这样你是不是可以想到设计另一种方法,满足对偶可行和对偶间隙为0,逐步改善原始不可行,下面就是对偶单纯形算法。


5. 对偶单纯形算法

  • 对偶单纯形是求解线性规划的另一个基本方法,它是根据对偶原理和单纯形法原理设计出来的,不要简单理解为是求对偶问题的单纯形法。
  • 基本思路:找出一个对偶问题的可行基,保持对偶问题为可行解的条件下,通过互补松弛性条件变换到原始问题的解看是否可行,可行就是最优解,不可行的时候保持对偶问题为可行解的情况下转移原始问题到另一个可行解。这部分不是本文的关键,所以有些简略,可自行参考资料配合上一篇单纯形的文章去理解,下面我写出来几个需要注意的点。
  • 初始表中一定要满足对偶问题可行,即检验数满足最优判别准则。
  • 对偶单纯形法与普通单纯形法的换基顺序不一样,普通单纯形法是先确定进基变量后确定出基变量,对偶单纯形法是先确定出基变量后确定进基变量。
  • 确定出基变量,取 b_l=\min \ \{b_i|b_i<0\}其对应变量 x_l 做为换出基的变量。
  • 确定入基变量:普通单纯形法的最小比值是 \min \ \{ \frac{\alpha_i}{\beta_{ik}}| \ \alpha_{ik}>0\}其目的是保证下一个原问题的基本解可行,对偶单纯形法的最小比值是 \min \limits_j \ \{\left| \frac{\lambda_j}{a_{lj}} \right||a_{lj}<0\} ,其目的是保证下一个对偶问题的基本解可行。


6. 原始-对偶算法

前面墨迹了这么久,终于到了本文的关键点了。原始-对偶算法的设计思路也是保证对偶可行,满足互补松弛性条件,逐步改善原始不可行性。

算法的基本过程如下:

下面我再用文字补充一下算法过程说明,顺便加上一些我自己的理解及细节。

(1)求(D)的一个初始可行解 y^0 :

  • c\ge0 时,取 y=0 即可;
  • c\ge0 不成立时,在原始问题(LP)中引进变量 x_{n+1} 和增加一个约束: x_1+x_2+\cdots+x_n+x_{n+1}=b_{m+1} 。同时取 c_{n+1}=0(LP^*) \min \ c^Tx+0*x_{n+1}\\ s.t. \ \begin{cases} A_i^Tx=b_i, i=1,2,\cdots,m\\ x_1+x_2+\cdots+x_n+x_{n+1}=b_{m+1}\\ x_i\ge0(i=1,2,\cdots,n,n+1) \end{cases} 可以看出这个问题和原问题具有相同的最优解,该问题的对偶问题为 (D^*)\max \ b^Ty+b_{m+1}y_{m+1}\\ s.t. \ \begin{cases} A_j^Ty+y_{m+1}\leq c_j,(j=1,2,\cdots,m)\\ y_{m+1} \leq0 \end{cases} (D^*) 有一个可行解 \begin{cases} y_i=0,(i=1,2,\cdots,m)\\ y_{m+1}=\min \ \{c_j|1\leq j \leq m \}<0 \end{cases}在这个寻找初始可行解的过程同样可以见到引入人工变量和原始对偶问题如何转换的情况,所以这两点trick一定要记住。

(2)求 y^k 的允许列集合 J_k=\{j:A_j^Ty^k=c_j\} x_i(c_i-A_i^Ty)=0这些列是满足最优解的互补松弛性条件的,后面还会结合迭代更新公式再次理解允许列集合。

(3)通过(LP)问题构造限定的原问题(RP)

如果原问题可行是满足下面的条件的:Ax=b\Leftrightarrow\sum_{j=1}^{n}{a_{ij}x_j}=b_i(i=1,2,\cdots,m)\Leftrightarrow\sum_{j\in J}{a_{ij}x_j}+\sum_{j\notin J}{a_{ij}x_j}=b_i\Leftrightarrow \sum_{j\in J}{a_{ij}x_j}=b_i

当时上课时老师提出一个问题为什么构造(RP)加一个非负项就可成立,我的理解是因为构造的问题是求 \min ,所以必须 x_i^a \ge0 ,如果允许 x_i^a<0 则构造的(RP)问题无界,此(RP)问题等价于原始问题(LP)。保证对偶可行,改善原始不可行,就是改善这个(RP)问题的最小值为0。

(4)把(RP)写成对偶形式(DRP),如果写参照前面的标准型往里套,把 \begin{bmatrix} x\\ x^a \end{bmatrix} 看成变量,目标函数添加 0*x ,有心的朋友不妨自己推推看,也是一个锻炼的过程。

(5)又要用到单纯形法求解了。可以得到(RP)问题的最优值 \xi_{opt} ,最优解 x^k 以及(DRP)的最优解 \bar y ,可能有人要问怎么得到(DRP)的最优解 \bar y ,这就要用到前面提到的互补松弛性条件的应用了,所以前面的铺垫都不是废话,还是很有用处的,一定要仔细看一遍。

  • 如果(RP)的最优值 \xi_{opt}=0 ,得到最优解,算法终止,这一点很好理解
  • 重点我来解释一下 \xi_{opt}>0 的情况,同时帮助理解一下允许列集合。先假定迭代公式 y^{k+1}=y^k+\theta \bar y ,则有 A^Ty^{k+1}=A^Ty^k+\theta A^T \bar yA^T \bar y=\begin{bmatrix} A_j^T \bar y\ \ (j \in J) ,\ A_j^T \bar y \ \ (j \notin J)\end{bmatrix}(RP)问题中的最优基里每个允许列,在下次迭代是还是保持允许的。 如果 A_j 是在(RP)的最优基里面,那么它的检验数为0,参考第四节写的公式先把(RP)问题写成 \begin{bmatrix} A^*,I \end{bmatrix} \begin{bmatrix} x\\x^a \end{bmatrix} =bz = \begin{bmatrix} 0,1 \end{bmatrix} \begin{bmatrix} x\\x^a\end{bmatrix} 的形式,方便推导。可得检验数 \begin{bmatrix}A^{*T}\\I \end{bmatrix} \bar y- \begin{bmatrix}0\\1 \end{bmatrix}= \begin{bmatrix} 0\\\lambda\end{bmatrix}A^*A 的允许列集合,所以有 A_j^T \bar y=0(j\in J) 满足(DRP)问题的可行域。 \begin{cases} A_j^Ty^{k+1}=A_j^Ty^k+\theta A_j^T \bar y = c_j \ \ \ j \in J\\ A_j^Ty^{k+1}=A_j^Ty^k+\theta A_j^T \bar y \ \ \ j \notin J \end{cases}如果所有的 A_j^T \bar y \leq0,j \notin J ,则 A_j^Ty^{k+1}<c_j ,说明允许列和非允许列没有变化,则下次迭代的(RP)问题还是相同的解,陷入死循环,而且 b^Ty^{k+1}=b^Ty^k+\theta b^T \bar y=b^Ty^k+ \theta \xi_{opt}>b^Ty^k ,所以问题无上界,原始问题不可行,算法不能进行下去。 如果 \exists j \notin J 使得 A_j^T \bar y>0 ,要维持 y^{k+1}=y^k+\theta \bar y 的可行性,需要使 A_j^Ty^{k+1}=A_j^Ty^k+\theta A_j^T \bar y \leq c_j 。即 \theta = \min \limits_{A_j^T\bar y>0,j\notin J}\{\frac{c_j-A_j^Ty^k}{A_j^T \bar y}\}

7. 结语

写到这里差不多该交代的都已经交代清楚了,这次没有写代码,我也没想好怎么编写这个算法的代码呢,结合到具体实际问题应用再考虑吧,关于原始对偶算法的应用有很多,比如应用到最短路问题就相当于反向的Dijkstra算法,这个也是可以推导出来的,我觉得这两篇文章已经够用了,受之予渔。





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


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


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

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

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

编辑于 2018-06-08

文章被以下专栏收录