單純形算法詳細解析

線性規划（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）的套路，即給定可行域的一個頂點，求值，沿一條邊到達下一個頂點，看是否能改善解，直到達到停止要求。

這里就要問幾個問題了：為什么最優解一定在頂點處？怎么得到頂點？怎么實現沿着一條邊到下一個頂點？什么時候停止？接下來，我們將一一解答這些問題，當解答完這些問題，單純形法也就顯而易見了。

1、凸多邊形最優解在頂點

考慮最小化問題，目標函數\(\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} \]

因此，總有一個頂點，他的目標函數值不比內部點差。

2、怎樣得到一個凸多邊形的頂點？

下面要說明的就是這樣一個定理：對於可行域方程組\(\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/