單純形算法詳細解析

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

這裏就要問幾個問題了：爲何最優解必定在頂點處？怎麼獲得頂點？怎麼實現沿着一條邊到下一個頂點？何時中止？接下來，咱們將一一解答這些問題，當解答完這些問題，單純形法也就顯而易見了。

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

考慮最小化問題，目標函數\(\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} \]
所以，總有一個頂點，他的目標函數值不比內部點差。

二、怎樣獲得一個凸多邊形的頂點？

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