Linear Inequality
System of Linear Inequalities
「一次不等式組」。許多道一次不等式同時成立。
{ 1 x + 2 y - 5 z < 5 { 2 x + 4 y + 6 z ≥ 1 { 3 x + 1 y + 7 z < 4
幾何意義是半空間交集,所有解形成一個凸多胞形,也可能形成開放區間、退化、空集合。請見本站文件「Half-plane」。
目前沒有演算法。大家更傾向討論Linear Programming。
有人適度乘上負號,調整成Ax > b的格式,套用高斯消去法,但是不知道正不正確。
Linear Programming
Mathematical Programming
「數學規劃」即是約束最佳化。已經發展成一系列問題:
linear programming min cᵀx subject to Ax ≤ b quadratic programming min 1/2 xᵀQx + cᵀx subject to Ax ≤ b second order cone programming min fᵀx subject to ‖Aᵢx + bᵢ‖ ≤ cᵢᵀx + dᵢ, Fx = g semidefinite programming min C●X subject to Aᵢ●X ≤ bᵢ, X ≽ 0 min tr(CX) subject to tr(AᵢX) ≤ bᵢ, X ≽ 0
LP ⊂ QP ⊂ SOCP ⊂ SDP
現代社會經常使用數學規劃,例如經濟交易、交通運輸、工業生產。數學規劃儘管不是電資科系的學習重點,然而卻是商管科系的必修專業。
經典數學領域亦可使用數學規劃,例如組合最佳化。數學規劃比起經典組合算法,有過之而無不及,計算時間更短。
Minimum Ratio Spanning Tree: Dinkelbach's Algorithm http://sunmoon-template.blogspot.com/2016/12/0-1-fractional-programming-dinkelbach.html
知名工具Gurobi、CPLEX、GAMS。大家把現實世界問題改寫成數學規劃問題,運用工具快速得到答案。
Linear Programming
「一次規劃」、「線性規劃」。一次函數的約束最佳化。目標函數、約束函數都是一次函數。
minimize x₁ + 2x₂ - x₃ + 1 subject to 3x₁ + 2x₂ + x₃ + 1 < 3 -4x₁ - 5x₂ + 6x₃ + 1 ≥ 4 7x₁ + 8x₂ - 9x₃ - 2 = 5
Linear Programming標準式
調整式子,美化式子。數學家的最愛。
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ ≤ 2 4x₁ + 5x₂ - 6x₃ ≤ -3 7x₁ + 8x₂ - 9x₃ ≤ 7 -7x₁ - 8x₂ + 9x₃ ≤ -7 x₁, x₂, x₃ ≥ 0
步驟很簡單,只是有點多。
一、目標函數,一律移除常數。 二、目標函數,一律求最大值。(式子視情況乘上-1。) 三、約束條件,一律將常數移項至右側。 四、約束條件,一律包含等於。 五、約束條件,一律是不等式。(一道等式拆成兩道不等式。) 六、約束條件,一律是小於等於。(式子視情況乘上-1。) 七、約束條件,一律是非負變數。(一個無界變數拆成兩個非負變數相減。)
一、目標函數,一律移除常數。
minimize x₁ + 2x₂ - x₃ subject to 3x₁ + 2x₂ + x₃ + 1 < 3 -4x₁ - 5x₂ + 6x₃ + 1 ≥ 4 7x₁ + 8x₂ - 9x₃ - 2 = 5
二、目標函數,一律求最大值。(式子視情況乘上-1。)
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ + 1 < 3 -4x₁ - 5x₂ + 6x₃ + 1 ≥ 4 7x₁ + 8x₂ - 9x₃ - 2 = 5
三、約束條件,一律將常數移項至右側。
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ < 2 -4x₁ - 5x₂ + 6x₃ ≥ 3 7x₁ + 8x₂ - 9x₃ = 7
四、約束條件,一律包含等於。
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ ≤ 2 -4x₁ - 5x₂ + 6x₃ ≥ 3 7x₁ + 8x₂ - 9x₃ = 7
五、約束條件,一律是不等式。(一道等式拆成兩道不等式。)
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ ≤ 2 -4x₁ - 5x₂ + 6x₃ ≥ 3 7x₁ + 8x₂ - 9x₃ ≤ 7 7x₁ + 8x₂ - 9x₃ ≥ 7
六、約束條件,一律是小於等於。(式子視情況乘上-1。)
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ ≤ 2 4x₁ + 5x₂ - 6x₃ ≤ -3 7x₁ + 8x₂ - 9x₃ ≤ 7 -7x₁ - 8x₂ + 9x₃ ≤ -7
七、約束條件,一律是非負變數。(一個無界變數拆成兩個非負變數相減。)
maximize -(x₁₊ - x₁₋) - 2(x₂₊ - x₂₋) + (x₃₊ - x₃₋) subject to 3(x₁₊ - x₁₋) + 2(x₂₊ - x₂₋) + (x₃₊ - x₃₋) ≤ 2 4(x₁₊ - x₁₋) + 5(x₂₊ - x₂₋) - 6(x₃₊ - x₃₋) ≤ -3 7(x₁₊ - x₁₋) + 8(x₂₊ - x₂₋) - 9(x₃₊ - x₃₋) ≤ 7 -7(x₁₊ - x₁₋) - 8(x₂₊ - x₂₋) + 9(x₃₊ - x₃₋) ≤ -7 x₁₊, x₁₋, x₂₊, x₂₋, x₃₊, x₃₋ ≥ 0
實際應用當中,變數通常就是非負,或者最佳解位置通常就是非負。這種時候,直接追加非負條件,不影響答案。大可不必拆變數。
maximize -x₁ - 2x₂ + x₃ subject to 3x₁ + 2x₂ + x₃ ≤ 2 4x₁ + 5x₂ - 6x₃ ≤ -3 7x₁ + 8x₂ - 9x₃ ≤ 7 -7x₁ - 8x₂ + 9x₃ ≤ -7 x₁, x₂, x₃ ≥ 0
至此得到標準式。
Linear Programming矩陣運算形式
標準式可以寫成矩陣運算。
max cᵀx subject to Ax ≤ b, x ≥ 0 [ -1 ] [ x₁ ] [ 3 2 1 ] [ 2 ] c = [ -2 ] x = [ x₂ ] A = [ 4 5 -6 ] b = [ -3 ] [ 1 ] [ x₃ ] [ 7 8 -9 ] [ 7 ] [ -7 -8 9 ] [ -7 ]
Linear Programming幾何意義
先備知識請見本站文件「Half-plane」。
可行解:滿足約束條件的解。 最佳解:滿足約束條件與目標函數的解。
(以三維空間為例) 變數:座標軸。 目標函數:平面,朝著法向量移動。 約束函數:半空間。 可行解:半空間交集。 最佳解:半空間交集的邊界,可能包含點、邊、面。
「Fundamental Theorem of Linear Programming」:可行解形成凸多胞形。最佳解位於頂點,可能連成一片。
因此得以使用貪心法。隨便一個可行解,逐步修改可行解,朝著目標函數的方向走,維持在界內,必定可以走到其中一個最佳解。
「Upper Bound Theorem for Polytope」:N個變數(維度)、M道不等式(刻面),可行解至多C(M-⌈N/2⌉, ⌊N/2⌋) + C(M-⌊N/2⌋-1, ⌈N/2⌉-1)個頂點。
「d-Step Conjecture」:N個變數(維度)、M道不等式(刻面),兩頂點的最短路徑至多M-N條邊。仍是懸案。
因此貪心法的時間複雜度仍是懸案。只知道至多是指數時間,不知道是否為多項式時間。
Projective Scaling Algorithm(投影縮放法)
行走於可行解內部。理論上是多項式時間,實務上極慢。
請見專著《Operations Research: Applications and Algorithms》的Karmarkar's Algorithm。
Simplex Algorithm(單形法)
行走於可行解的頂點與邊。理論上是指數時間,實務上極快。
Linear Programming: Simplex Algorithm
演算法
單形法有許多種解讀方式。
《Operations Research: Applications and Algorithms》仿效高斯消去法 《Introduction to Algorithms》仿效鬆弛法 《Introduction to Linear Optimization》仿效最陡下降法
此處採用高斯消去法的觀點。因為我只熟悉這個觀點。
一、一次規劃化作一次方程組。 二、先找可行解:藉由消去法,逐步消滅負變數。 三、再找最佳解:藉由消去法,逐步拉升最大值。
標準式(Standard Form)
為了容易實作程式碼,首先調整成標準式。
maximize x₁ + 2x₂ subject to 8x₁ + 3x₂ ≤ 24 2x₁ - 9x₂ ≤ -18 3x₁ + 4x₂ ≤ 12 x₁, x₂ ≥ 0
等式(Slack Variables)
目標函數:最大值視作新變數。
約束條件:添加非負變數。小於等於被補足成等於。
x₁ + 2x₂ = z 8x₁ + 3x₂ + s₁ = 24 2x₁ - 9x₂ + s₂ = -18 3x₁ + 4x₂ + s₃ = 12 x₁, x₂, s₁, s₂, s₃ ≥ 0
消去法(Pivoting)
橫條的四則運算,不影響可行解、最佳解。
以任一位置為基準,實施消去法,不影響可行解、最佳解。
x₁ + 2x₂ = z 8x₁+ 3x₂+ s₁ = 24 2x₁ - 9x₂ + s₂ = -18 3x₁ + 4x₂ + s₃ = 12 x₁, x₂, s₁, s₂, s₃ ≥ 0 -13/3 x₁ - 2/3 s₁ = z-16 8/3 x₁+ x₂+ 1/3 s₁ = 8 78/3 x₁ + 9/3 s₁ + s₂ = 54 -23/3 x₁ - 4/3 s₁ + s₃ = -20 x₁, x₂, s₁, s₂, s₃ ≥ 0
設定初始值(Basic Solution)
原始變數設為0,添加變數設為b,最大值設為0。
實施消去法,一個變數從b變0,一個變數從0變b。
nonbasic basic maximum nonbasic basic maximum variable variable value variable variable value x₁ = 0 s₁ = 24 0 pivoting x₁ = 0 x₂ = 8 16 x₂ = 0 s₂ = -18 -------> s₁ = 0 s₂ = 54 s₃ = 12 s₃ = -20
尋找可行解(Basic Feasible Solution)
教科書是大M法、兩階段法,實務上是消去法。
變數不可以是負數,必須通通調整成非負數。
一、常數項b(末端直條),尋找負數。 回、沒負數,即是可行解。立即結束。 回、有負數,需要調整成非負數,繼續下一步。 二、該橫條,尋找負數。 回、沒負數,即是無解。立即結束。 回、有負數,那還有救。以該變數為準,實施消去法。 該變數、常數項b,調整成非負數了。
↓ x₁ + 2x₂ = z 8x₁ + 3x₂ + s₁ = 24 2x₁ - 9x₂ + s₂ =-183x₁ + 4x₂ + s₃ = 12 x₁, x₂, s₁, s₂, s₃ ≥ 0 ↓ x₁ + 2x₂ = z 8x₁ + 3x₂ + s₁ = 24 2x₁- 9x₂+ s₂ = -18 ← 3x₁ + 4x₂ + s₃ = 12 x₁, x₂, s₁, s₂, s₃ ≥ 0 13/9 x₁ + 2/9 s₂ = z-4 26/3 x₁ + s₁ + 1/3 s₂ = 18 -2/9 x₁ + x₂ - 1/9 s₂ = 2 35/9 x₁ + 4/9 s₂ + s₃ = 4 x₁, x₂, s₁, s₂, s₃ ≥ 0
尋找最佳解(Basic Optimal Solution)
以消去法逐步拉升最大值。
一、目標函數(首端橫條),尋找正數。 回、沒正數,即是最佳解。立即結束。 回、有正數,繼續下一步,以便拉升最大值。 二、該直條,尋找正數。 回、沒正數,即是最佳解無限大。立即結束。 回、有正數,可拉升最大值。以該變數為準,實施消去法。
13/9 x₁+ 2/9 s₂ = z-4 ← 26/3 x₁ + s₁ + 1/3 s₂ = 18 -2/9 x₁ + x₂ - 1/9 s₂ = 2 35/9 x₁ + 4/9 s₂ + s₃ = 4 x₁, x₂, s₁, s₂, s₃ ≥ 0 ↓ 13/9 x₁ + 2/9 s₂ = z-4 ← 26/3 x₁ + s₁ + 1/3 s₂ = 18 -2/9 x₁ + x₂ - 1/9 s₂ = 235/9 x₁+ 4/9 s₂ + s₃ = 4 x₁, x₂, s₁, s₂, s₃ ≥ 0 + 2/35 s₂ - 13/35 s₃ = z-192/35 s₁ - 23/35 s₂ - 78/35 s₃ = 318/35 x₂ - 3/35 s₂ + 2/35 s₃ = 78/35 x₁ + 4/35 s₂ + 9/35 s₃ = 36/35 x₁, x₂, s₁, s₂, s₃ ≥ 0
+ 2/35 s₂- 13/35 s₃= z-192/35 ← s₁ - 23/35 s₂ - 78/35 s₃ = 318/35 x₂ - 3/35 s₂ + 2/35 s₃ = 78/35 x₁ + 4/35 s₂ + 9/35 s₃ = 36/35 x₁, x₂, s₁, s₂, s₃ ≥ 0 ↓ + 2/35 s₂ - 13/35 s₃ = z-192/35 s₁ - 23/35 s₂ - 78/35 s₃ = 318/35 x₂ - 3/35 s₂ + 2/35 s₃ = 78/35 x₁ + 4/35 s₂+ 9/35 s₃= 36/35 x₁, x₂, s₁, s₂, s₃ ≥ 0 -1/2 x₁ + - 1/2 s₂ = z-6 23/4 x₁ + + s₁ - 3/4 s₂ = 15 3/4 x₁ + x₂ + 1/4 s₂ = 3 35/4 x₁ + + 9/4 s₂ = 9 x₁, x₂, s₁, s₂, s₃ ≥ 0 optimum = 6
鬼打牆(Cycling)
實施消去法,極少數情況下,最大值持續不變,不斷鬼打牆。
計算順序(Pivot Rule)
請見《Pivoting rules for the revised simplex algorithm》。
精挑細選消去法標的,避免鬼打牆。方法很多,列舉幾個:
一、大小順序(Dantzig's rule): (可能鬼打牆) 首端橫條,從中選擇係數最大者。 該直條,從中選擇拉升最多者。 二、索引順序(Bland's rule): (時間略增) 首端橫條,從中選擇索引值最小者。 該直條,從中選擇拉升最多者、然後索引值最小者。 三、字典順序(lexicographic rule): (時間略增) 首端橫條,從中任選。 所有橫條各自縮放,將該直條係數變成一。 該直條,從中選擇橫條字典順序最小者。 四、次數順序(Zadeh's rule): (空間略增) 首端橫條,從中選擇使用次數最少者。 該直條,從中選擇拉升最多者。 五、貪心順序(greatest increment rule): (時間暴增) 選擇拉升最多者。
程式碼
此處採用大小順序,可能鬼打牆。
為了方便實作迴圈,目標函數從首端挪至末端。
為了節省時間空間,不紀錄單元向量。再建立兩個陣列,紀錄變數位置。
[ 8 3 24 ] [ 26/3 1/3 18 ] [ 2 -9 -18 ] pivoting [ -2/9 -1/9 2 ] [ 3 4 12 ] -------> [ 35/9 4/9 4 ] [ 1 2 0 ] [ 13/9 2/9 -4 ] (objective) x₁ x₂ b x₁ s₂ b
Integer Linear Programming
Integer Linear Programming
「整數一次規劃」。所有變數必須是整數。整數解。
整數一次規劃有個特別之處:限制變多了,解變少了,求解卻變難了!一次規劃是P問題,整數一次規劃卻是NP問題。
演算法請見專著《Operations Research: Applications and Algorithms》。
使用方式請見專著《Integer Linear Programming in Computational and Systems Biology: An Entry-Level Text and Course》。
Total Unimodular Matrix and Integer Matrix => P problem https://math.stackexchange.com/questions/64461/ integer-matrices-with-integer-inverses => Unimodular matrix https://math.stackexchange.com/questions/19528/
Branch and Bound(分支定界法)
整數一次規劃:每次將其中一個變數xᵢ的「邊界(上界或下界)」確定下來。因為x的大小與f(x)的大小有直接關聯,所以藉由調整每個變數xᵢ的上下界,夾擠出f(x)的極值。
整數規劃:每次將其中一個變數xᵢ的「確切數值」確定下來。因為x的大小與f(x)的大小沒有直接關聯,所以採用了比較複雜的策略:一開始設定f(x)的極值的下界是零,隨著已知變數越來越多,逐步增加f(x)的極值的下界,以逼近f(x)的極值。至於f(x)的極值的上界,並沒有用來尋找極值,主要是用來阻止冗餘的分支。
Cutting Plane Method(切割平面法)
整數一次規劃:分數解所在的約束條件,任取其中一個,擷取小數部分,使之小於1(小於等於0),形成新的約束條件。追加此約束條件,雖然影響分數解,但是不影響整數解。追加約束條件、修正分數解,不斷重複直到形成整數解。
整數規劃:梯度下降法沿著梯度方向走一步,而切割平面法則是答案範圍之內隨便挑一點(例如重心),沿著梯度垂直方向割一刀,縮小答案範圍。
Quadratic Programming🚧
Quadratic Programming
「二次規劃」。目標函數是二次函數,約束函數是一次函數。
NP-hard。大家轉為討論正定矩陣。
大家沿用內點法。
Newton–Lagrange Method(牛頓拉格朗日乘數法)
拉格朗日乘數,牛頓法。
請見專著《Numerical Optimization》SQP章節。
min f(x) subject to c(x) = 0 L(x,λ) = f(x) + λᵀc(x) solve { d/dx L(x,λ) = 0 { d/dλ L(x,λ) = 0 solve { f′(x) + c′(x)ᵀλ = 0 { c(x) = 0 solve [ f′ c′ᵀ] [ x ] = [ 0 ] [ c 0 ] [ λ ] [ 0 ] newton method: solve [ f″ c′ᵀ] [ dx ] = [ f′+c′ᵀ ] [ c′ 0 ] [ dλ ] [ c ] let [ x_next ] = [ x ] - [ dx ] [ λ_next ] [ λ ] [ dλ ]
Uzawa's Algorithm(宇澤法)
拉格朗日乘數,共軛梯度法。
https://en.wikipedia.org/wiki/Uzawa_iteration
min 1/2 xᵀAx - bᵀx subject to Cx = d L(x,λ) = (1/2 xᵀAx - bᵀx) + λᵀ(Cx - d) solve { d/dx L(x,λ) = 0 { d/dλ L(x,λ) = 0 solve { Ax + Cᵀλ - b = 0 { Cx - d = 0 solve [ A Cᵀ] [ x ] = [ b ] [ C 0 ] [ λ ] [ d ]
Active Set Method(有效集合法)
https://blog.csdn.net/dymodi/article/details/50358278
http://www.cs.nthu.edu.tw/~cherung/teaching/2011cs5321/handout8.pdf
找到一個可行解。挑選幾個限制條件。KKT,求解。取消一個限制,KKT,求解。兩解相減作為行進方向,步伐大小盡量大,更換限制條件。不斷行進,直到乘數非負。
二次函數,偏微分等於零,即是一次方程組,可以用高斯消去法求解。
Sequential Quadratic Programming(循序二次規劃)
約束最佳化化作二次規劃。利用拉格朗日乘數、KKT條件。
Wolfe's Modified Simplex Method(單形法)
二次規劃化作一次規劃。
拉格朗日乘數、單形法、限制某些變數不可同時非零。
請見專著《Operations Research: Applications and Algorithms》。
Trust Region Method(信賴區域法)
最佳化化作二次最佳化。
當前地點走一步,其二階泰勒近似,在步伐大小半徑範圍之內(信賴區域),求得極值位置,作為下一處。二次函數的極值位置必定位於邊界之上,因此步伐大小其實不會短少。
根據爬升幅度,決定要不要走,幅度夠大才走。根據爬升幅度,調整步伐大小,幅度越大則下次步伐越大。
若極值很難求,可以改用梯度、曲率倒數的加權平均數,作為前進方向。宛如擬牛頓法。