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。大家轉為討論正定矩陣。
大家沿用內點法。
equality-constrained quadratic programming
「約束條件都是等式的二次規劃」。約束函數是等式。
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(宇澤法)
拉格朗日乘數,共軛梯度法。
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 ]
trust region(信賴區域)
最佳化化作二次最佳化。
地面勘查類型的最佳化演算法,可以分成兩類。
一、直線搜尋(撞名了):人為指定步伐方向暨步伐大小,例如一次微分(不動點遞推法)(梯度下降法)、二次微分(牛頓法)、二次微分改成其他類似的東西(擬牛頓法)、混搭(Levenberg–Marquardt algorithm)。
二、信賴區域:自動尋找步伐方向暨步伐大小,成為最佳化問題。每走一步,就解一個最佳化問題。
直線搜尋可以視作信賴區域的特例。每走一步,就解一個最佳化問題,其答案是公式解:一次微分、二次微分、二次微分改成其他類似的東西。
言歸正傳。自動尋找步伐方向暨步伐大小的手法如下:
一、當前地點,函數求得一階泰勒近似。
二、在指定半徑範圍之內(信賴區域),求得極值位置,作為下一個地點。
三、根據爬升幅度,決定要不要走,幅度夠大才走。
四、根據爬升幅度,調整半徑範圍(間接調整步伐大小),幅度越大則下次半徑範圍越大(步伐大小越大)。
sequential quadratic programming(循序二次規劃)
約束最佳化化作二次規劃。
約束最佳化,套用拉格朗日乘數、KKT條件,形成多變數函數最佳化。演算法採用牛頓法,步伐方向暨步伐大小是一次微分(梯度)除以二次微分。
多變數函數二次微分得到矩陣(Hessian matrix)。由於零很多,通常沒有反矩陣,只好改成其他類似的的東西。宛如擬牛頓法。
改為自動尋找步伐方向暨步伐大小。步伐必須滿足約束條件。仿照信賴區間,採用一階泰勒展開,下一個地點必須在約束條件之內。
目標函數的零次項是常數,可以省略。一次項可以從拉格朗日函數改為原本函數。當約束條件是等式,根據新約束條件,一次項當中的原約束函數部分形成常數,可以省略。當約束條件是不等式,透過添加非負變數slack variable硬是改成等式。但是非負變數要怎麼處理呢?沒講!到頭來還是只能處理等式約束條件?
active set method(有效集合法)
二次規劃化作約束條件都是等式的二次規劃。
找到一個可行解。挑選幾個限制條件。KKT,求解。取消一個限制,KKT,求解。兩解相減作為行進方向,步伐大小盡量大,更換限制條件。不斷行進,直到乘數非負。
二次函數,偏微分等於零,即是一次方程組,可以用高斯消去法求解。
Wolfe's modified simplex method(單形法.改)
二次規劃化作約束最佳化。
拉格朗日乘數、單形法、限制某些變數不可同時非零。
請見專著《Operations Research: Applications and Algorithms》。