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🚧

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》。