[ simplex algorithm ] 線性規劃 - 單純形法

線性規劃簡介

線性規劃是在線性約束條件下尋找某個目標線性函數的最大、最小值，例如：
$$
\begin{array}{lrl}
    \min& -3x_1+2x_2 & \\
    \mathrm{subject~to}
        &4x_1-6x_2&\le 5\\
        &5x_1-2x_2&\ge 7\\
        &x_1&\ge 0.
    \end{array}
$$

在最佳化問題中，線性規劃是非常重要的領域，很多的演算法問題都可以轉換成線性規劃問題來解決，例如：最大網路流、最大匹配等。因此研究如何系統化的解決線性規劃問題是非常重要的事情

線性規劃標準形

所有的線性規劃問題都可以轉換成下列標準型：

$$
\begin{array}{lrl}
    \max &\sum_{j=1}^n c_j \times x_j &\\
    \mathrm{subject~to}
        &\sum_{j=1}^n A_{i,j} \times x_j &\le b_j ,\; \forall i=1\sim m\\
        &x_j& \geq 0, \; \forall  j = 1\sim n
    \end{array}
$$

一般來說大多數解線性規劃的演算法都會要求將問題轉換成標準形，因此務必要理解其轉換方法。轉換成標準形的步驟如下：

1. 將最小化改成最大化：
   只要改變目標函數的正負號，即可將最小化問題改成標準型設定的最大化問題，也就是說$$\min \sum_{j=1}^n c_j \times x_j$$等價於$$\max \sum_{j=1}^n -c_j \times x_j$$
2. 去除等式：
   如果約束條件中存在等式，將其轉化為兩個分別為大於等於及小於等於的不等式取代，例如：
   $$x_1+x_2=5$$等價於$$x_1+x_2 \le 5\\x_1+x_2 \ge 5$$
3. 去除大於等於：
   如果約束條件中存在大於等於約束，將約束兩邊取負即可，例如：
   $$x_1+x_2 \ge 5$$ 等價於 $$-x_1-x_2 \le -5$$
4. 去除自由變數：
   如果未知數 $x_i$ 沒有非負約束，我們說 $x_i$ 是自由變數。加入新變量 $x_i'$，並用 $x_i-x_i'$替換原來的變量 $x_i$，並增加$x_i,x_i' \ge 0$的條件即可

標準化範例：

$$
\begin{array}{lrl}
    \min& -x_1+2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&= 5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
首先將最小化改成最大化：
$$
\begin{array}{lrl}
    \max& x_1-2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&= 5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
接著去除等式：
$$
\begin{array}{lrl}
    \max& x_1-2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&\le 5\\
        &x_1+x_2&\ge 5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
去除大於等於：
$$
\begin{array}{lrl}
    \max& x_1-2x_2 & \\
    \mathrm{subject~to}
        &x_1+x_2&\le 5\\
        &-x_1-x_2&\le -5\\
        &x_1-x_2&\le 2\\
        &x_2&\ge 0.
    \end{array}
$$
最後去除自由變數，將$x_1$用$x_1-x_3$取代：
$$
\begin{array}{lrl}
    \max& x_1-2x_2-x_3 & \\
    \mathrm{subject~to}
        &x_1+x_2-x_3&\le 5\\
        &-x_1-x_2+x_3&\le -5\\
        &x_1-x_2-x_3&\le 2\\
        &x_1,x_2,x_3&\ge 0.
    \end{array}
$$

單純形法

輸入input

本人實作的單純形法其輸入是一個被稱為單純形表的 $m+1 \times n+1$ 矩陣$A$，標準形中每個項目的值都可以對應到矩陣$A$上：
$$
\begin{array}{lrl}
    \max &\sum_{j=1}^n A_{0,j} \times x_j &\\
    \mathrm{subject~to}
        &\sum_{j=1}^n A_{i,j} \times x_j &\le A_{i,0} ,\; \forall i=1\sim m\\
        &x_j& \geq 0, \; \forall  j = 1\sim n
    \end{array}
$$
可以發現：

1. $A_{0,1} \sim A_{0,n}$對應標準形中的$c_1 \sim c_n$
2. $A_{i,0}$對應標準形中的$b_i ,\; \forall i=1\sim m$

為了方便計算請設$A_{0,0}=0$

上面標準化範例寫成單純形表為以下形式：

$$
\begin{bmatrix}
0& 1 &-2 &-1 \\
5 &1& 1 &-1 \\
-5 &-1 &-1 &1 \\
2 &1 &-1& -1
\end{bmatrix}
$$

輸出output

輸出為一個$size=n+1$的vector $ans$，其中：

1. $ans[0]$為線性規劃的解，也就是最大值
2. $ans[1] \sim ans[n]$為變數$x_1 \sim x_n$的其中一組滿足條件的答案
3. 如果無解或是解無限大則$ans$的$size=0$
4. 注意如果答案是 $0$ 由於浮點數誤差的關係可能會產生 $``-0"$ 這種東西，要小心處理

上面標準化範例的解如下：

$$\{0.5,\; 3.5,\; 1.5,\; 0\}$$

意即在 $x_1=3.5,\;x_2=1.5,\;x_3=0$ 的情況下會有最大值 $x_1-2x_2-x_3=0.5$

程式碼

使用C++11

原理

有空時會補上，先讓我好好休息吧！終於可以正常睡覺了

[ Implement Queue using Stacks ] 用堆疊實作佇列

一般在C++ STL中queue都是用deque實作的，而deque用和vector類似的倍增法來實作(稍微多了一些操作導致常數有點大但是每次操作時間較為平均)，雖然夠用但有些時候這不是一個好選擇。
大二資料結構期中考時，考了一題叫你用stack實作queue的方法，我那時候想了一下想了一個超帥的作法，需要用到兩個stack，我把它成為stack A和B。

• push:
  push時直接把資料push進B中，$\ord{1}$
• pop:
  這是個問題，我們把A的top當成是queue的front，所以直接把A pop就可以了，如果A是empty的話，就把B的資料依序pop然後push進A中，在進行該操作，均攤$\ord{1}$
• front:
  如pop所述，把A的top當成是queue的front，如果A是empty再另做處理，$\ord{1}$
• back:
  和front相反，把B的top當成是queue的back，如果B是empty再另做處理，$\ord{1}$

以下為實作過程，以c++ vector代替stack：

[ Median-of-Medians Algorithm ] 中位數演算法

這是一個可以在保證線性時間(c++ std::nth_element是隨機演算法)找出一個序列中第k大元素的演算法，網路上已經有不少教學，但是很多人都認為常數太大因此缺乏實作。

教學文在此：http://tmt514-blog.logdown.com/posts/484313-divide-and-conquer-method-iii-find-the-median

其實我高中時就想要試著去時做看看，但是因為那時的程式能力太差的關係，做出來的東西一直有bug，後來去忙其他事情後就被我忘掉了。最近因為學長面試有被問到一樣的問題跑來問我，才慢慢想起來有一份沒寫完的code，於是今天抱著不管怎樣都要寫出來的精神把他寫完了：

[ C++ std::sort python implement ] C++std::sort Python實作

最近有一個很靠北的課需要用python寫一些C++很容易做到的東西。
這次我們被要求寫一個quick sort，但是他要求的做法實在太糟糕了，於是我就參考了C++ std::sort來實做了一個quick sort(不包含heap sort的部分)：

[ inorder postorder construct tree ] 用中序、後序建樹保證$\ord{N}$的方法

昨天我學長因為要面試，所以努力地刷leetcode的題目，寫到了一題:

• 106. Construct Binary Tree from Inorder and Postorder Traversal  

他雖然AC了但是用的方法不太好，因此跑來問我有沒有看起來很帥速度又快的方法。

因為網路上的code都用一些很慢的方法來建樹，很多都$\ord{N^2}$，雖然也有看似$\ord{N}$的方法，但是用到像是unordered_map之類常數很大也不保證複雜度是$\ord{1}$(codeforces會卡內建hash)。

因為我查到的code都太糟糕了，因此我決定自己寫一個複雜度保證為$\ord{N}$又好寫的方法。

首先是找$root$的部分，因為有給後序的關係很容易就是到是誰了，但是找左右子樹就會出問題。經過觀察，只需要在dfs的過程用stack紀錄一些特別點就可以好好地維護左右子樹。

假設現在dfs在點$u$，stack紀錄的點就是從$root$到$u$的路徑所有點中，除了那些左小孩也在該路徑中的點之外的所有點，有點複雜看個圖會比較明白，紅色的點就是記錄在stack中的點。

至於為什麼記錄這些點就可以在dfs時判斷現在是不是NULL的節點，以及如果給的是preorder的情況就交給讀者思考了
以下附上程式碼:

[ Source code beautifier / syntax highlighter ] 在網頁/blog中插入彩色程式碼

先看看結果吧:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

#define x first #define y second #include<bits/stdc++.h> using namespace std; #define X(){\ sdfgsdfg;\ sdfgsdfg;\ } int main(){ //asdfasdfdfsdfghd\ asdfasdfasdf\ asdfasdfsdf wchar_t wc; cout<<"Jinkela"<<'\n'; cout<<R"jinkela( 7122)jinkela"; cout<<L"adsfasdf"<<endl; return 0; /* asdf */ }

這是使用http://hilite.me/ Style=monokai的結果，個人覺得效果不錯，只是raw string和一些比較難實作的東西沒有支援而已，其他的都還算可以

我已經把這個網址加到我的學習連結裡面了，C/C++ syntax highlighter (Style選monokai)那個。用法就是貼上程式碼，設定好語言和style，把產生的html貼在你想貼的位置，蠻簡單的

我在blog中如果程式碼量比較少我覺得沒必要加入模板也會用這個方法來貼code

[ Steiner tree problem in graphs ] 斯坦納樹

斯坦納樹問題是一個世界知名的NP-hard問題。在圖論上的斯坦納樹是對於一張無向圖$G=(V,E)$以及一個點集合$P \subseteq V$，通常會稱$P$集合為$terminal \; set$，對於每條邊$e=(u,v) \in E$，令$w(e)$表示它的權重。我們的目標是要從$G$的所有子圖中找出一棵生成樹$T=(V',E')$，使得$P \subseteq V'$且$\sum_{e \in E'} w(e)$最小。

簡單來說就是在圖$G$上找一棵子樹，可以把$P$中的點連通起來，且邊權總和最小

如果我們枚舉所有子圖，對每個子圖做最小生成樹演算法，就一定可以找到斯坦納樹，但是複雜度是$\ord{(\abs E + \abs V log \abs V ) \times 2^{\abs V}}$，非常糟糕。

如果$w(e)>0 ,e \in E$，且$\abs P \ll \abs V$，我們可以找到一個動態規劃的方法：
令$dp[S][i]$表示以點$i$為根，以$S \subseteq P$為$terminal \; set$構造出來的斯坦納樹，這樣我們最後的答案就會是$dp[P][u \in P]$

狀態轉移式可以寫成這樣

• $dp[S][i]=min(dp[T][j]+dp[S-T][j]+dis(i,j)\; : \; j \in V,T \subset S)$
  $dis(i,j)$表示$i \sim j$的最短路徑

 

任兩點間的最短路徑可以用floyd在$\ord {\abs V^3}$預先算出來，狀態有$2^{\abs P}\times \abs V$個，狀態轉移為$\ord{\abs V \times 枚舉子集合的時間}$，因此總複雜度為$\ord{\abs V^3+2^{\abs P} \times \abs V^2 \times 枚舉子集合的時間 }$

其中 $2^{\abs P} \times 枚舉子集合的時間$ 只是粗略的計算，實際上它是
$$\sum_{i=1}^{\abs P} \binom{\abs P}{i} \times (2^i -1) \simeq \sum_{i=0}^{\abs P} \binom{\abs P}{i} \times 2^i = (1+2)^{\abs P} = 3^{\abs P}$$因此總複雜度可以表示為$\ord{V^3+3^{\abs P} \times \abs V^2}$，但是這其實還可以優化，令$H[j] = min(dp[T][j]+dp[S-T][j] \; : \;T \subset S)$
則$dp[S][i]=min(H[j]+dis(i,j)\; : \;j \in \abs V)$
$H$是可以被預先算出來的，因此總複雜度就降為$\ord{\abs V^3 + \abs V 3^{\abs P}+\abs V^2 2^{\abs P}}$
以下附上程式碼:

有的時候圖是稀疏圖，也就是$\ord V=\ord E$，這種時候用這種算法效率其實不是很好，我們可以在dp的過程才用一些單源最短路徑算法算出最短路徑，這樣複雜度可以變成$$\ord{\abs V 3^{\abs P}+ShortestPath(G) 2^{\abs P}}$$其中$ShortestPath(G)$是在圖$G$中計算最短路徑的時間，用dijkstra的話是$\ord{\abs E+\abs V log \abs V}$，這裡我用SPFA實作: