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

知名工具GurobiCPLEXGAMS。大家把現實世界問題改寫成數學規劃問題,運用工具快速得到答案。

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₂      = -18
3x₁ + 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₂      =   2
35/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

時間複雜度

每步需時O(NM)。

最差步數不明。有一種很差的情況是「Klee–Minty Cube」,需要走2ᴺ - 1步。最差時間複雜度是指數時間。

平均步數是M步。平均時間複雜度是多項式時間。

現今已知的計算順序,時間複雜度如上所述;尚未發明的計算順序,也許有更好的時間複雜度。單形法的最佳計算順序仍是懸案,單形法的時間複雜度仍是懸案。

UVa 10498 ICPC 7584 7836

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(Under Construction!)

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(信賴區域法)

最佳化化作二次最佳化。

當前地點走一步,其二階泰勒近似,在步伐大小半徑範圍之內(信賴區域),求得極值位置,作為下一處。二次函數的極值位置必定位於邊界之上,因此步伐大小其實不會短少。

根據爬升幅度,決定要不要走,幅度夠大才走。根據爬升幅度,調整步伐大小,幅度越大則下次步伐越大。

若極值很難求,可以改用梯度、曲率倒數的加權平均數,作為前進方向。宛如擬牛頓法。