Linear Optimization

Linear Optimization

「線性最佳化」。線性函數最佳化。線性函數求極值。

輸出變數從一個變成了多個,衍生兩個問題:

大小:輸出如何比大小?定義向量長度:所有元素的平方和。

梯度:線性函數的梯度?定義矩陣微分:對各個輸入變數微分。

矩陣微分(矩陣梯度)

一個線性函數,對各個變數微分。得到多個導數。

                              ∂f         ∂f         ∂f     
f(x) = (x₀ + x₁ + x₂)         ――― = 1    ――― = 1    ――― = 1
                              ∂x₀        ∂x₁        ∂x₂    

                              ∂f         ∂f         ∂f      
f(x) = (c₀x₀ + c₁x₁ + c₂x₂)   ――― = c₀   ――― = c₁   ――― = c₂
                              ∂x₀        ∂x₁        ∂x₂     

一個線性函數,對一整個向量微分。將導數們從上到下依序併成一個向量。

                                  [ x₀ ]
f(x) = (c₀x₀ + c₁x₁ + c₂x₂)   x = [ x₁ ]
                                  [ x₂ ]
  ∂f   [ ∂f/∂x₀ ]   [ c₀ ]
  ―― = [ ∂f/∂x₁ ] = [ c₁ ]
  ∂x   [ ∂f/∂x₂ ]   [ c₂ ]
∂   T                                        d        
―― y x = y                             --->  ―― xy = y
∂x                                           dx       

∂   T                                        d        
―― x y = y                             --->  ―― yx = y
∂x                                           dx       

∂   T                                        d   2     
―― x x = 2x                            --->  ―― x  = 2x
∂x                                           dx        
                      when A is
∂   T            T    symmetric              d    2      
―― x A x = (A + A ) x ========= 2Ax    --->  ―― ax  = 2ax
∂x                                           dx          

多個線性函數,依序併成一個向量,對整個向量微分。將導數們從左到右依序併成一個矩陣。

       [ f₀(x) ]   [ a₀x₀ + a₁x₁ + a₂x₂ ]   [ a₀ a₁ a₂ ]       
       [ f₁(x) ]   [ b₀x₀ + b₁x₁ + b₂x₂ ]   [ b₀ b₁ b₂ ] [ x₀ ]
F(x) = [ f₂(x) ] = [ c₀x₀ + c₁x₁ + c₂x₂ ] = [ c₀ c₁ c₂ ] [ x₁ ]
       [ f₃(x) ]   [ d₀x₀ + d₁x₁ + d₂x₂ ]   [ d₀ d₁ d₂ ] [ x₂ ]
       [ f₄(x) ]   [ e₀x₀ + e₁x₁ + e₂x₂ ]   [ e₀ e₁ e₂ ]       
                                                 A         x   

  ∂F     ∂f₀ ∂f₁ ∂f₂ ∂f₃ ∂f₄     [ a₀ b₀ c₀ d₀ e₀ ]       
  ―― = [ ――― ――― ――― ――― ――― ] = [ a₁ b₁ c₁ d₁ e₁ ]
  ∂x     ∂x  ∂x  ∂x  ∂x  ∂x      [ a₂ b₂ c₂ d₂ e₂ ]       
              ∂/∂x A x                   Aᵀ               
∂         T                                  d        
―― A x = A                             --->  ―― ax = a
∂x                                           dx       

∂   T                                        d        
―― A x = A                             --->  ―― ax = a
∂x                                           dx       

∂               T                            d              
―― (A x + b) = A                       --->  ―― (ax + b) = a
∂x                                           dx             

∂   T      ∂    T         T  T   T           d          
―― y A x = ―― (y A) x = (y A) = A y    --->  ―― axy = ay
∂x         ∂x                                dx         

多個線性函數,視作一個線性函數,輸入輸出是多個變數。對其微分,即是矩陣微分。

向量長度(向量範數)

定義長度的平方。向量所有元素的平方總和。即是向量自己與自己的點積。

定義微分連鎖律。因為矩陣複合是從右到左,所以連鎖律設定成從右到左。

定義長度的平方

     T      def₁: squared norm     2
(A x) (A x) ================== ‖Ax‖ 

定義連鎖律

∂      2 def₂: chain rule ∂     ∂       2    T            T   
―― ‖Ax‖  ================ ―― Ax ――― (Ax) = (A )(2Ax) = 2 A A x
∂x                        ∂x    ∂Ax                           

這兩個定義的合理性,可由嚴謹的推導過程證明:
                                            T
∂       T        ∂   T T      ∂   T  T     A A 必定對稱    T   
―― (A x) (A x) = ―― x A A x = ―― x (A A) x ============ 2 A A x
∂x               ∂x           ∂x                               

當A退化為I,依然合理:

∂   T    ∂     2  
―― x x = ―― ‖x‖ = 2x
∂x       ∂x       

∂   T    ∂      T       ∂      2     T
―― x x = ―― (Ix) (Ix) = ―― ‖Ix‖ = 2 I I x = 2x
∂x       ∂x             ∂x

範例:解線性方程組Ax = b

                 2
minimize ‖Ax - b‖

∂          2
―― ‖Ax - b‖ = 0                     「一次微分等於零」的地方是極值、鞍點
∂x                                   二次函數、開口向上,必是最小值
  ∂
[ ―― (Ax - b) ] [ 2(Ax - b) ] = 0   微分連鎖律
  ∂x
 T
A  [ 2(Ax - b) ] = 0                微分

 T       T
A A x = A b                         同除以2、展開、移項

      T -1 T
x = (A A) A b                       Normal Equation
min f(x) = ‖Ax - b‖²   =>   ∇f(x) = AᵀAx - Aᵀb = 0

結果宛如向量投影公式:b投影到A,求得投影量!

平方誤差盡量小=垂直投影!打通了分析與幾何兩大領域!

亦得用於線性迴歸。

範例:解線性方程組Ax = 0

             2                 2
minimize ‖Ax‖  , subject to ‖x‖ = 1

             2         2
minimize ‖Ax‖  - λ (‖x‖ - 1)     Lagrange's multiplier

∂        2         2
―― [ ‖Ax‖  - λ (‖x‖ - 1) ] = 0   「一次微分等於零」的地方是極值、鞍點
∂x

   T
2 A A x - 2 λ x = 0              微分

 T                                       T
A A x = λ x                      正解是 A A 最小的特徵值

亦得用於直線擬合。

範例:Rayleigh Quotient

                  2
minimize ‖Ax - λx‖

∂           2
―― ‖Ax - λx‖ = 0                       「一次微分等於零」的地方是極值、鞍點
∂λ                                      二次函數、開口向上,必是最小值
  ∂
[ ―― (Ax - λx) ] ∙ [ 2(Ax - λx) ] = 0   微分連鎖律
  ∂λ

[ -x ] ∙ [ 2(Ax - λx) ] = 0             微分

 T       T
x A x = x λ x                           同除以-2、展開、移項

     xᵀ A x
λ = ――――――――                            Rayleigh Quotient
      xᵀ x

結果宛如向量投影公式:Ax投影到x,求得投影量!

範例:Rayleigh Quotient導數

∂     xᵀ A x
―― [ ―――――――― ] = 0  and  ‖x‖² = xᵀ x ≠ 0
∂x     xᵀ x

 (A + Aᵀ) x     (xᵀ A x) 2x
―――――――――――― - ――――――――――――― = 0    分式微分
    xᵀ x          (xᵀ x)²

                xᵀ A x 
(A + Aᵀ) x = 2 ―――――――― x           移項,同乘以 xᵀ x,重新整理
                 xᵀ x

(A + Aᵀ) x = 2 λ x                  虛擬特徵值

 A + Aᵀ
――――――― x = λ x                     移項,形成特徵向量的格式
   2

Semidefinite Optimization(Under Construction!)

Semidefinite Optimization

「半定最佳化」。雖然名稱不分正負、涵蓋了零,但是一般是指正定矩陣最佳化。

二次型:xᵀAx。
正定矩陣:xᵀAx > 0,不討論x = 0。
對稱正定矩陣:xᵀAx > 0,不討論x = 0。A恰好對稱。

正定矩陣仿照了ax² > 0的概念,幾何意義是位於X軸上方、開口朝上、碰觸原點的拋物線。

ax² > 0推廣到xᵀAx > 0之後,輸入變數有許多個。例如兩個變數的時候,就是仿照ax₁² + bx₂² + cx₁x₂ + dx₂x₁ > 0的概念,幾何意義是橢圓拋物面,等高線是同心橢圓。

【註:輸入變數超過兩個的時候,橢圓軸線仍然互相垂直?我不太確定。如果你知道麻煩告訴我。】

正定矩陣是橢圓拋物面函數,極值位於橢圓軸線上方。按道理應該有公式解,可是我沒見過。正定矩陣亦是凸函數,得以套用高速的最佳化演算法,例如梯度下降法。

多個橢圓拋物面函數相加、乘上正倍率,仍是橢圓拋物面函數。

[X1,X2] = meshgrid(-10:1:10,-10:1:10); Y = (X1 .* X1 .* 10) + (X2 .* X2 .* 6) + (X1 .* X2 .* 2); surf(X1,X2,Y)

ax² > 0推廣到xᵀAx > 0之後,迸出了「對稱」的概念。

對稱正定矩陣,比正定矩陣的性質更強。特徵向量即是橢圓軸線,特徵值平方根倒數即是橢圓軸長。

對稱正定矩陣,極值擁有公式解(共軛梯度法?)。最大(小)值位於最大(小)的特徵值的特徵向量上方。

[X1,X2] = meshgrid(-10:1:10,-10:1:10); Y = (X1 .* X1 .* 10) + (X2 .* X2 .* 6) + (X1 .* X2 .* 8); surf(X1,X2,Y)
[E,D] = eig([2 1; 1 3]) [E,D] = eig([2 0; 2 3])
[X1,X2] = meshgrid(-10:1:10,-10:1:10); Y = (X1 .* X1 .* 2) + (X2 .* X2 .* 3) + (X1 .* X2 .* 2); surf(X1,X2,Y)
ContourPlot[(x * x * 2) + (y * y * 3) + (x * y * 2), {x, -4, 4}, {y, -4, 4}] Plot3D[(x * x * 2) + (y * y * 3) + (x * y * 2), {x, -4, 4}, {y, -4, 4}]
Plot3D[(x * x * 2) + (y * y * 3) + (x * y * 2), {x, -6, 6}, {y, -6, 6}, Boxed -> False, Axes -> False, MeshFunctions -> {#3 &}, Mesh -> 15, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &) ]

Positive-definite = Convex【尚待確認】

正定函數明明就是凸函數,卻有人稱作單調函數。哭夭啦。

https://math.stackexchange.com/questions/2550058/

function f is monotone iff (y-x) dot (f(y)-f(x)) >= 0
                                     ^^^^^^^^^^^
                                     f(y-x) when f is linear/affine

linear/affine function is monotone iff positive-semidefinite

vector field is monotone iff jacobian matrix is positive-semidefinite everywhe
re.

範例:二次型的極值

對稱正定矩陣,最大(小)值位於最大(小)的特徵值的特徵向量上方。利用正定最佳化來證明。

          T                    2
minimize x A x , subject to ‖x‖ = 1  (A is symmetric positive-definite)

          T             2
minimize x A x - λ ( ‖x‖ - 1 )      Lagrange multiplier

∂     T             2
―― [ x A x - λ ( ‖x‖ - 1 ) ] = 0   「一次微分等於零」的地方是極值、鞍點
∂x                                  二次型,必為極值

2 A x - 2 λ x = 0                   展開

A x = λ x                           移項,形成特徵向量的格式

以上結論是針對最小值。針對最大值,結果亦類似。

Multivariate Optimization(Under Construction!)

Multivariate Optimization

「多變量最佳化」。輸入與輸出有許多個變數。

n->n求根   有人說是eigenvalue
           變成companion matrix? Cayley-Hamilton Theorem?
http://www.fing.edu.uy/if/cursos/fiscomp/extras/numrec/book/f9.pdf

n->n牛頓法 根據n維版本泰勒展開式而得 每個維度不可分開處理 需要jacobian
            物理意義不明
http://fourier.eng.hmc.edu/e176/lectures/NM/node21.html

n->n割線法  通常被視作類牛頓法  一樣有jacobian
https://en.wikipedia.org/wiki/Broyden's_method
n->n最佳化 到底是什麼東西? 矩陣多項式?
https://en.wikipedia.org/wiki/Vector_optimization

n->n梯度下降法 物理意義完全不明 最佳化對象到底是啥? 但是也有 jacobian
https://en.wikipedia.org/wiki/Gradient_descent#Solution_of_a_non-linear_system

n->n牛頓法解特徵值
http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter3.pdf

多函數的一次微分

多個函數,依序拼成一個向量。多個函數對向量微分,將導數們從左到右依序拼成一個矩陣。

    [ f₀ ]
F = [ f₁ ]
    [ f₂ ]

∂F     ∂f₀ ∂f₁ ∂f₂
―― = [ ――― ――― ――― ] = [ ∂f₀/∂x  ∂f₁/∂x  ∂f₂/∂x  ]
∂x     ∂x  ∂x  ∂x
                       [ ∂f₀/∂x₀ ∂f₁/∂x₀ ∂f₂/∂x₀ ]
                     = [ ∂f₀/∂x₁ ∂f₁/∂x₁ ∂f₂/∂x₁ ]
                       [ ∂f₀/∂x₂ ∂f₁/∂x₂ ∂f₂/∂x₂ ]

                     = J(F)

該矩陣是「多函數的一階導數矩陣Jacobian Matrix」。此處的定義仿照了線性函數微分,但是數學家的定義卻是此矩陣的轉置。應該是古人沒想清楚?

單函數的二次微分

一次微分(函數的梯度):函數對向量微分,得到向量。

二次微分(函數的梯度的梯度):再次對向量微分,得到矩陣。

     ∂f   [ ∂f/∂x₀ ]
∇f = ―― = [ ∂f/∂x₁ ]
     ∂x   [ ∂f/∂x₂ ]

 2    ∂∇f     ∂(∂f/∂x₀)  ∂(∂f/∂x₁)  ∂(∂f/∂x₂)  
∇ f = ――― = [ ―――――――――  ―――――――――  ――――――――― ]
       ∂x         ∂x         ∂x         ∂x     

            [ ∂²f/∂x₀∂x₀ ∂²f/∂x₁∂x₀ ∂²f/∂x₂∂x₀ ]
          = [ ∂²f/∂x₀∂x₁ ∂²f/∂x₁∂x₁ ∂²f/∂x₂∂x₁ ]
            [ ∂²f/∂x₀∂x₂ ∂²f/∂x₁∂x₂ ∂²f/∂x₂∂x₂ ]

          = H(f) = J(∇f)

該矩陣是「單函數的二階導數矩陣Hessian matrix」。

該矩陣恰巧是「多函數的一階導數矩陣Jacobian Matrix」代入「函數的梯度」,然而這個觀點似乎沒有什麼用處。

演算法(Newton Method)

牛頓法求根:從座標(x, f(x))的地方畫切線,與座標軸相交之處,做為下一個x。

∇f(x) (x - xnext) = f(x)

牛頓法多函數求根:多個函數形成一個向量。所有函數一齊實施牛頓法求根。
             
Jᵀ(f(x)) (x - xnext) = f(x)

牛頓法求極值:極值(與鞍點)出現在一次微分等於零的地方。
       預先微分,再套用牛頓法多函數求根。

Hᵀ(f(x)) (x - xnext) = ∇f(x)

牛頓法多函數求極值:像這種要求,我這輩子沒見過。

Taylor Series

                    1               1          2
f(x + Δx) = f(x) +  ―― f'(x) Δx  +  —— f"(x) Δx  + ...
                    1!              2!

                    1               1
f(x + Δx) = f(x) +  ——   Jᵀ  Δx  +  —— Δxᵀ Hᵀ Δx + ...
                    1!              2!
  n
Δx projects onto basis, which is function gradient.
textbook use different definition. textbook use transpose.

我的Jacobian Matrix有別於數學家的版本,因此我的牛頓法、泰勒展開式也有別於數學家的版本。我不知道我的版本是否合理。也許泰勒展開式應該解讀成投影?

Linear Programming(Under Construction!)

前言

最佳化當中,「規劃Programming」是指約束最佳化,而且約束條件只有不等式。等式可以通通換成不等式:以大於等於、小於等於的兩道不等式,來取代一道等式。

已經發展出許多類型:

LP ⊂ QP ⊂ SOCP ⊂ SDP
linear programming
min cᵀx   s.t. Ax ≤ b, x ≥ 0

quadratic programming
min 1/2 xᵀQx + cᵀx   s.t. Ax ≤ b

second order cone programming
min fᵀx   s.t. ‖Aᵢx + bᵢ‖ ≤ cᵢᵀx + dᵢ, Fx = g

semidefinite programming
min C●X   s.t. Aᵢ●X ≤ bᵢ, X ≽ 0 (positive semi-definite)
min tr(CX)   s.t. tr(AᵢX) ≤ bᵢ, X ≽ 0

Linear Programming

「線性規劃」。目標函數、約束條件都是線性函數,只考慮第一象限。

幾何意義,請參考計算幾何「Half-plane Intersection」。目標函數:一個方向向量。約束條件:半平面。可行解:半平面交集。最佳解:半平面交集的頂點、邊、面。

因為解是凸函數,所以可以設計極快的演算法!

https://reference.wolfram.com/language/ref/RegionPlot3D.html
max f(x)     ---> min -f(x)
max min f(x) ---> max g(x) s.t. f(x) > g(x)

現代社會經常使用線性規劃。舉凡經濟交易、交通運輸、工業生產,都能看到線性規劃的應用。雖然線性規劃不是電資科系的學習重點,卻是商管科系的專業必修。

一些經典數學領域,例如組合最佳化、排程理論,除了傳統的組合算法,也能用線性規劃解題,有過之而無不及。

專著如《Understanding and Using Linear Programming》。工具如LINDO

Simplex Algorithm(單形法)

http://ocw.mit.edu/courses/aeronautics-and-astronautics/16-410-principles-of-autonomy-and-decision-making-fall-2010/lecture-notes/MIT16_410F10_lec17.pdf

一、變成線性方程組。目標函數,預設答案,視作等式。約束條件,不等式添加變數,成為等式。

二、求可行解。取原點做為可行解:原始變數設為0,添加變數設為b。如果原點不是可行解:添加變數為負數,有兩種解法。

二甲、兩階段法。添加變數為負值者,替其添加暫時變數。奢望暫時變數是0,故新增目標函數:最小化暫時變數加總。以高斯消去法消去暫時變數,使之為0,令新目標函數變成約束條件。最後刪去暫時變數。

二乙、大M法。添加變數為負值者,替其添加暫時變數。目標函數,減去暫時變數,並且乘上巨大係數,使得最佳解的暫時變數必為零。

三、求最佳解。貪心法,藉由高斯消去法,等價調整約束條件,逐步提高目標函數值。幾何意義是:可行解,沿邊走,朝向目標函數的方向。

三甲、高維度的情況下,運氣非常不好時,可能走進一大片鞍點,在山腰平原上鬼打牆。解法是小小擾動b,摧毀鞍點。

【待補程式碼】

每步需時O(N(M+N))。

N個變數(維度)、M道不等式(刻面),可行解至多C(M,N)個頂點。根據「Upper Bound Theorem」,可行解至多M⌊N/2⌋個頂點。

單形法至多走M⌊N/2⌋步,最差時間複雜度是指數時間。

單形法平均走M步,平均時間複雜度是多項式時間。

一種很差的情況是「Klee-Minty Cube」,需要走2ᴺ - 1步,但是遭遇機率極低。

UVa 10498 ICPC 7584 7836

Cutting Plane Method(切割平面法)

答案範圍內部隨便挑一點(例如重心),求得梯度,沿梯度垂直方向割一刀,縮小答案範圍。

梯度下降法則是沿著梯度方向走一步。

The center of gravity method就是挑重心。

凸函數的時間複雜度O(N log 1/ϵ)。

Interior Algorithm(內點法)
Barrier Method(屏障法)

梯度下降法。約束條件取log,以Regularization添入原式,效果是擠壓邊界、調整路徑。

f(x) + log(g(x))

Semidefinite Programming

Semidefinite Programming

「半定規劃」。請參考專著《Numerical Optimization》

Convex Programming

Convex Programming

「凸規劃」。目標函數與約束函數是凸函數。

上述章節通通都是凸規劃的特例。

請參考專著《Convex Optimization》

判斷一個函數是否為凸函數,是NP-complete問題。於是有人改弦易轍,排列組合常見運算符號、常見函數,判斷是否為凸函數,甚至判斷要用哪個演算法。

請參考專論《Disciplined Convex Programming》

Alternating Direction Method of Multipliers

http://stanford.edu/class/ee364b/lectures/admm_slides.pdf

目標函數是凸函數們的線性函數,約束條件是線性函數。

min u f(x) + v g(x)   subject to a x + b y = 0

Integer Programming

Integer Programming

「整數規劃」。解必須是整數。

Integer Programming                     整數規劃
Integer Linear Programming              整數線性規劃
Integer Semidefinite Programming        整數半定規劃
Mixed Integer Programming               混和整數規劃(有些變數不必是整數)
Mixed Integer Linear Programming        混和整數線性規劃
Mixed Integer Semidefinite Programming  混和整數半定規劃

工具如GurobiCPLEX

整數規劃有個特別之處:限制變多了,解變少了,求解卻變難了!線性規劃是P問題,整數線性規劃卻是NP問題。原因不明。

Branch and Bound(分支定界法)

這是一個過時的演算法,觀念陳舊,已被淘汰。

n個變數x1…xn,找f(x)的可行解。

在線性規劃問題當中,b&b每次將其中一個變數xi的「邊界(可能是上限或者是下限)」給確定下來。由於x的大小與f(x)的大小有直接關連,所以藉由調整每個變數xi的上下限,就能夾擠出f(x)的極值。

https://www.youtube.com/watch?v=BbrZsG7zesE

在離散組合問題當中,b&b又不一樣了。b&b每次將其中一個變數xi的「確切數值」給確定下來。由於x的大小與f(x)的大小沒有直接關連,所以採用了比較複雜的策略:一開始設定f(x)的極值的下限是零,隨著已知變數越來越多,逐步增加f(x)的極值的下限,以逼近f(x)的極值。至於f(x)的極值的上限,並沒有用來尋找極值,主要是用來阻止多餘的分支。

https://www.youtube.com/watch?v=nN4K8xA8ShM

在離散組合問題當中,b&b跟state space search基本相同。唯一的差異在於,作業研究(工管)講到b&b的時候,通常會簡化矩陣的row和column,只拿不為零的數值來分支。人工智慧(資工)講到state space search的時候,是使用原本矩陣。

Convex Duality(Under Construction!)

Lagrange Duality

可微函數的對偶。製造新維度,窮舉斜率λ,求最高截距。

http://www.eng.newcastle.edu.au/eecs/cdsc/books/cce/Slides/Duality.pdf
http://www.onmyphd.com/?p=duality.theory
http://www.cnblogs.com/90zeng/p/Lagrange_duality.html

Legendre Duality(Convex Conjugate)

凸函數的對偶。畫個切線求截距。

動態規劃有一個技巧是:最佳解位於直線們的包絡線,經過點線對偶,最佳解位於凸包切線的截距。正是此對偶。

maximum: largest value
minimum: smallest value
supremum: least upperbound function
infimum: most lowerbound function
epigraph: upper area (of convex function)
hypograph: below area (of concave function)
https://en.wikipedia.org/wiki/Fenchel's_duality_theorem
let f is convex, g is concave, f*(y) = sup { xᵀy - f(x) }
inf { f(x) - g(x) } = sup { g*(y) - f*(y) }
D f(x) = [D f*(y)]⁻¹ 
sup { ⟪x,∇f(x)⟫ - f(x) }
Fenchel inequality: f(x) + f*(y) ≥ ⟪x,y⟫
Moreau decomposition: x = prox f(x) + prox f*(x)
infimal convolution: (f ◇ g)(x) = inf { f(y) + g(x − y) }
Fenchel duality: (f ◇ g)* = f* + g* , (f + g)* = f* ◇ g*.

Linear Programming Duality

線性函數的對偶。線性規劃必有拉格朗日對偶。

https://people.ok.ubc.ca/bauschke/Research/68.pdf
max cx   subject to Ax ≤ b
min by   subject to Aᵀy ≥ c