單純形演算法的原理和示例實現
作者簡介： @磊磊 本科河北工業大學電子科學與技術專業，現在中科院半導體研究所讀碩士，目前研一，方向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; beta = BN; 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];endend
推薦閱讀：