Linear Optimization

Linear Optimization

「線性最佳化」。一次函數最佳化。

此處討論稍微複雜的版本,輸入輸出是向量。

f(x) = ax + b推廣成F(x) = Ax + b。

Plot[2 x + 1, {x, -3, 3}, PlotRange -> {-3,3}, AspectRatio -> 1] VectorPlot[{2x+y+1, -x+y}, {x, -3, 3}, {y, -3, 3}]

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

一、如何比大小?各個輸出變數的平方和‖F(x)‖²。

二、如何求微分?分別對各個輸入變數微分F′(x)。

矩陣微分(矩陣梯度)

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

                              ∂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₃ ]

        ∂         [ ∂/∂x₁ f(x) ]   [ c₁ ]
f′(x) = ―― f(x) = [ ∂/∂x₂ f(x) ] = [ c₂ ]
        ∂x        [ ∂/∂x₃ f(x) ]   [ c₃ ]
∂                                            d        
―― yᵀx = y                             --->  ―― xy = y
∂x                                           dx       

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

∂                                            d         
―― xᵀx = 2x                            --->  ―― x² = 2x
∂x                                           dx        
                     when A is
∂                    symmetric               d           
―― xᵀAx = (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′(x) = ―― F(x) = [ f₁′(x) f₂′(x) f₃′(x) f₄′(x) f₅′(x) ]
        ∂x        [   |      |      |      |      |    ]

                  [ a₁ b₁ c₁ d₁ e₁ ]
                = [ a₂ b₂ c₂ d₂ e₂ ]
                  [ a₃ b₃ c₃ d₃ e₃ ]
                          Aᵀ               
∂                                            d        
―― A x = Aᵀ                            --->  ―― ax = a
∂x                                           dx       

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

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

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

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

向量長度(向量範數)

定義向量長度平方。向量所有元素的平方和。向量自身點積。

‖x‖² := xᵀx       = x₁² + x₂² + x₃² + ... + xₙ²
                 
‖Ax‖² = (Ax)ᵀ(Ax) = (Ax₁)² + (Ax₂)² + (Ax₃)² + ... + (Axₙ)²

向量長度的微分

定義微分連鎖律。矩陣連乘是從右到左,連鎖律也是從右到左。

∂             ∂       ∂      
―― F(G(x)) := ―― G(x) ―――― F(G(x))
∂x            ∂x      ∂G(x)

∂          ∂    ∂                                     
―― ‖x‖²  = ―― x ―― (x)² = (1)(2x) = 2x
∂x         ∂x   ∂x                              

∂          ∂     ∂                                     
―― ‖Ax‖² = ―― Ax ――― (Ax)² = (Aᵀ)(2Ax) = 2AᵀAx
∂x         ∂x    ∂Ax                              

定義的合理性,可由嚴謹的推導過程證明:

                                                  AᵀA is always      
∂          ∂              ∂           ∂             symmetric        
―― ‖Ax‖² = ―― (Ax)ᵀ(Ax) = ―― xᵀAᵀAx = ―― xᵀ(AᵀA)x ============= 2AᵀAx
∂x         ∂x             ∂x          ∂x                             

極值唯一性

‖Ax + b‖²有唯一極值,除非遇到特例。

‖⋅‖²向量長度平方是拋物面函數,也是凸函數,唯一極值,沒有鞍點。「一次微分等於零」的地方必是極值、不是鞍點。得以手工推導公式解,不需要特殊演算法。

g = Plot3D[Norm[{2x+y+1, -x+y}]^2, {x, -3, 3}, {y, -3, 3}, PlotRange -> {-3,3} , BoxRatios->{1,1,1}, ClippingStyle -> None, Mesh -> None, AxesOrigin->{0,0,0}, Boxed -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]
g = Plot3D[Norm[{2x+y+1, -x+y}]^2, {x, -3, 3}, {y, -3, 3}, MeshStyle -> LightGray, AxesOrigin->{0,0,0}, Boxed -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]; img = ImageResize[Rasterize[g, "Image", ImageResolution -> 72*3], Scaled[1/3]]; ImageCrop[RemoveBackground[img, {{{0, 0}}, 0.1}]]

數學公式(Normal Equation)

minimize ‖Ax + b‖²

∂
―― ‖Ax + b‖² = 0                    「一次微分等於零」的地方是極值、鞍點
∂x                                   二次函數、開口向上,必是最小值

  ∂
[ ―― (Ax + b) ] [ 2(Ax + b) ] = 0   微分連鎖律
  ∂x

Aᵀ [ 2(Ax + b) ] = 0                微分

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

x = -(Aᵀ A)⁻¹ Aᵀ b                  Normal Equation加上負號

範例:線性迴歸

minimize ‖Ax - b‖²

∂
―― ‖Ax - b‖² = 0                    「一次微分等於零」的地方是極值、鞍點
∂x                                   二次函數、開口向上,必是最小值

  ∂
[ ―― (Ax - b) ] [ 2(Ax - b) ] = 0   微分連鎖律
  ∂x

Aᵀ [ 2(Ax - b) ] = 0                微分

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

x = (Aᵀ A)⁻¹ Aᵀ b                   Normal Equation

範例:直線擬合

minimize ‖Ax‖²   subject to ‖x‖² = 1

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

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

2 AᵀA x - 2 λ x = 0                微分
                                
AᵀA x = λ x                        正解是 AᵀA 最小的特徵值

單位向量變換之後盡量短(長)=特徵向量!

Quadratic Optimization(Under Construction!)

Quadratic Optimization

「二次最佳化」。二次函數最佳化。

此處討論稍微複雜的版本,輸入是向量、輸出是數值。

f(x) = ax² + bx + c推廣成f(x) = xᵀAx + bᵀx + c。

Plot[3*x*x + -2*x + 1, {x, -3, 3}, PlotRange -> {-3,3}, AspectRatio -> 1] Plot3D[(3x-2y)*x + (x+y)*y + x + 1, {x, -3, 3}, {y, -3, 3}, PlotRange -> {-3, 3}, BoxRatios->{1,1,1}, ClippingStyle -> None, Mesh -> None, AxesOrigin->{0,0,0}, Boxed -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]

Quadratic Form

ax²是經過原點的拋物線。

xᵀAx是經過原點的拋物面。稱作「二次型」。

g = Plot3D[-x*x+y*y+x*y, {x, -3, 3}, {y, -3, 3}, PlotRange -> {-3, 3} , BoxRatios->{1,1,1}, ClippingStyle -> None, Mesh -> None, Axes->None, Boxed -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &), RegionFunction -> Function[{x, y, z}, x*x + y*y + z*z <= 9]]; h = ParametricPlot3D[{{u, 0, 0}, {0, u, 0}, {0, 0, u}}, {u, -3, 3}, PlotStyle -> Table[{Black, Thickness[0.01]}, 3]]; img = ImageResize[Rasterize[Show[g,h], "Image", ImageResolution -> 72*3], Scaled[1/3]]; ImageCrop[RemoveBackground[img, {{{0, 0}}, 0.05}]]

Quadratic Form: Positive Definite Matrix

當a > 0,則ax² > 0。除了x = 0是例外。

ax² > 0是位於X軸上方、開口朝上、碰觸原點的拋物線。

當A是正定矩陣,則xᵀAx > 0。除了x = 0是例外。

xᵀAx > 0是位於上半部、開口朝上、碰觸原點的橢圓拋物面。

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

等高線是同心橢圓。

Plot3D[x*x*3+y*y-x*y, {x, -3, 3}, {y, -3, 3}, Axes->None, Boxed -> False, MeshStyle->{Red,Thick}, MeshFunctions -> (#3 &), Mesh -> {Table[{j}, {j, 0, 50, 2}]}, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &) ]

高峻處、低窪處位於橢圓軸線上方。

【註:各種尺寸的同心圓的極值,連成一線。我不知道這種數學概念的正式名稱。似乎有許多神奇的應用方式。】

ContourPlot[x*x*3+y*y-x*y, {x, -3, 3}, {y, -3, 3}, Frame -> False, ContourShading -> None, Contours -> 21]
[E,D] = eig([3 -2; 1 1]) [E,D] = eig([3 -0.5; -0.5 1])

Quadratic Form: Symmetric Positive Definite Matrix

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

矩陣的「對稱」並不等於函數圖形的「對稱」。矩陣的「對稱」導致特徵向量互相垂直。白話解釋是「內心端正」。

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

高峻處、低窪處,位於最大、最小特徵值的特徵向量上方。

minimize xᵀAx   subject to ‖x‖² = 1   A is symmetric positive-definite

minimize xᵀAx - λ ( ‖x‖² - 1 )        Lagrange multiplier

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

2 A x - 2 λ x = 0                     微分

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

正定矩陣A,計算(A+Aᵀ)/2,得到對稱正定矩陣(A+Aᵀ)/2。兩者的二次型函數圖形一模一樣!因此教科書只討論對稱正定矩陣。

極值唯一性

f(x) = xᵀAx + bᵀx + c是否有唯一極值,取決於二次項xᵀAx。

二次項xᵀAx是曲面、一次項bᵀx是平面、零次項c是常數。

Plot3D[{x*x*3+y*y-x*y, x, 1}, {x, -3, 3}, {y, -3, 3}, PlotRange->{-3,3}, BoxRatios->{1,1,1}, ClippingStyle -> None, Mesh->None, Axes->None, Boxed -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]

二次項xᵀAx是否有唯一極值,取決於矩陣A。

正定矩陣有唯一最小值、負定矩陣有唯一最大值。

演算法(矩陣求解)

一次微分等於零的地方是極值或鞍點,而正定矩陣沒有鞍點。f′(x) = (A + Aᵀ)x + b = 0,問題化作矩陣求解。

高斯消去法O(N³)。Jacobi Method O(N²K)。N是維度。

演算法(Steepest Descent Method)

「最陡下降法」。給定地點x,給定前進方向d,正向逆向皆可,請一步走到該方向的最低處xnext

畫成圖形就是等高線的切線d和切點xnext

ContourPlot[(3x-2y)*x + (x+y)*y + x + 1, {x, -3, 3}, {y, -3, 3}, Frame -> False, ContourShading -> None, Contours -> 21]

正定矩陣A化做對稱正定矩陣(A+Aᵀ)/2,二次函數圖形一模一樣。方便起見,以下直接討論對稱正定矩陣A。f′(x) = 2Ax + b。

推得步伐大小-0.5 dot(f′(x), d) / dot(Ad, d)。

current position is x
current direction is d
find next position xnext such that F(xnext) is minimum along the path.

{ xnext = x + d ⋅ step
{ dot( f′(xnext) , d ) = 0

dot( f′(x + d ⋅ step)     , d ) = 0
dot( 2A(x + d ⋅ step) + b , d ) = 0
dot( 2Ax + 2Ad ⋅ step + b , d ) = 0
dot( 2Ax + b + 2Ad ⋅ step , d ) = 0
dot( f′(x)   + 2Ad ⋅ step , d ) = 0
dot(f′(x), d) + dot(2Ad, d) ⋅ step = 0
step = -dot(f′(x), d) / dot(2Ad, d)
step = -0.5 dot(f′(x), d) / dot(Ad, d)

梯度下降法(當前點的法線方向)、最陡下降法(下一點的切線方向)、座標下降法(座標軸方向)。全部介紹完畢。

演算法(Conjugate Direction Method)

「共軛方向法」。承上,仔細規定每一步的前進方向。

一、梯度反方向(朝向低處)前後垂直dot(dᵢ, dᵢ₊₁) = 0:逐步趨近最小值。
二、互相垂直dot(dᵢ, dⱼ) = 0:最多N步,無法抵達最小值。
三、互相傍A垂直dot(Adᵢ, dⱼ) = 0:最多N步,必定抵達最小值。

「共軛方向法」是第三種方式。正確性證明:

一、最小值位置x*的性質:位於一次微分等於零的地方。
f′(x*) = 2Ax* + b = 0

二、誤差e = x* - x的性質:誤差(變換後)恰是梯度一半。
Ae = Ax* - Ax = (Ax* + 0.5b) - (Ax + 0.5b) = -0.5 f′(x)

三、前進方向d的性質:垂直於梯度、誤差(變換後)。
dot(dᵢ, Aeᵢ₊₁) = -0.5 dot(dᵢ, f′(xᵢ₊₁)) = 0

四、初始誤差e₀ = x* - x₀的性質:步伐大小a,加總N步。
e₀ = a₀d₀ + a₁d₁ + ... + aN-1dN-1

五、前進方向互相傍A垂直,推得步伐大小無誤。
e₀ = a₀d₀ + a₁d₁ + ... + aN-1dN-1
Ae₀ = A(a₀d₀ + a₁d₁ + ...)
dot(d₀, Ae₀) = dot(d₀, A(a₀d₀ + a₁d₁ + ...))
dot(d₀, Ae₀) = a₀ dot(d₀, Ad₀) + a₁ dot(d₀, Ad₁) + ...
dot(d₀, Ae₀) = a₀ dot(d₀, Ad₀)
a₀ = dot(d₀, Ae₀) / dot(d₀, Ad₀)
a₀ = -0.5 dot(d₀, f′(x₀)) / dot(d₀, Ad₀)
a₀ = step₀
similarly, a₁ = step₁ , a₂ = step₂, ...

註:四和五理應從步伐大小開始推導,式子順序上下顛倒。
  由於那樣太複雜了,只好改成驗算的方式。

挑選N個前進方向,一種方式是Krylov Subspace聯合Gram-Schmidt Orthonormalization,將r, Ar, A²r, A³r, ..., Aᴺ⁻¹r扳成互相傍A垂直。r可以是任意值。

共軛方向法和高斯消去法的時間複雜度都是O(N³),但是共軛方向法可以逐步趨近最小值。實務上,共軛方向法只取足夠精確度,得以取消最後幾步,進而縮短時間,優於高斯消去法。

共軛方向法跟共軛毫無關聯!論文作者誤以為「傍A垂直A-orthogonal」可稱做「共軛conjugate」。

演算法(Conjugate Gradient Method)

承上,挑選N個前進方向,採用每一處的梯度反方向(朝向低處),將-f′(x)們扳成互相傍A垂直,邊走邊扳。

想深入鑽研細節的讀者,可以參考《An Introduction to the Conjugate Gradient Method Without the Agonizing Pain》

演算法詳細步驟,借用一下維基百科的圖片:

r := -f′(x)。p := d。α := step。

上圖當中,f(x) = 1/2 xᵀAx - bᵀx,f′(x) = Ax - b,目的是解Ax = b。如果要最佳化f(x) = xᵀAx + bᵀx + c,記得先將A乘以2、將b變號,再套用上圖。

Quadratic Form: Adjacency Matrix

圖論「Adjacency Matrix」,推廣成加權版本。

xᵀWx是兩兩相乘、乘上權重、通通加總!權重是W的元素。

一點:wx²。兩點:w₁₁x₁² + w₂₂x₂² + w₁₂x₁x₂ + w₂₁x₂x₁。三點:w₁₁x₁² + w₂₂x₂² + w₃₃x₃² + w₁₂x₁x₂ + w₂₁x₂x₁ + w₁₃x₁x₃ + w₃₁x₃x₁ + w₂₃x₂x₃ + w₃₂x₃x₂。

應用廣泛。

1. Multivariate Normal Distribution。Covariance Matrix。
2. Principal Component Analysis。Representation Matrix。

Quadratic Form: Laplacian Matrix

L = D - W  (diagonal of D is sum of W)
x*Lx = sum Wij (xi - xj)²

圖論「Laplacian Matrix」,推廣成加權版本。

xᵀLx兩兩相減平方、乘上權重、通通加總!權重是W的元素。

一點:0。兩點:w₁₂(x₁-x₂)²。三點:w₁₂(x₁-x₂)² + w₁₃(x₁-x₃)² + w₂₃(x₂-x₃)²。

W必須是對稱矩陣,才能構成相減平方。W不是對稱矩陣,那就改成對稱矩陣,Wsym = (W + Wᵀ) / 2或Wsym = D-1/2 W D-1/2

W有很多種選擇、很多種功能。

12_Parameterization1.pdf
smoothing (averaging): Wij = 1
distance preserving: Wij = |xi - xj|???
harmonic (angle preserving): Wij = 0.5 (cot(a) + cot(b)) ab是邊ij對頂兩角
mean-value: Wij = (tan(a/2) + tan(b/2)) / 2|xi - xj|    ab是邊ij與點j旁邊兩角

高峻處、低窪處,位於最大、最小特徵值的特徵向量上方。

L is symmetric positive semidefinite.
orthonormal eigenvectors, non-negative eigenvalues
min solution is eigenvalue zero, eigenvector [1,1,1,1,1,1,...]

特徵函數,然並卵。得到互相垂直的分量,進而做為座標軸。

curve:    t->x(t) t->y(t) t=0~1    two function
polyline: i->x[i] i->y[i] i=0~n-1  two array
surface:  uv->x(u,v) uv->y(u,v) uv->z(u,v)  three 2d function
mesh:  uv merge to i

eigenfunction of laplace operator of curve:    t->f(i)   infinite functions
                                                         (plot as nodal sets)
eigenvector   of laplace operator of polyline: i->f(i)   n vectors
                                                         (sin/cos wave)
eigenfunction/vector plot on curve/polyline as various colors. [-1,+1]->[0,255]
nodal set of eigenfunction -----> meaningless?
nodal set of Fiedler vector ----> look likes rings

應用廣泛。

1. Energy Minimization。請參考本篇文章。
2. Embedding。請參考本站文件「Embedding」。
3. Mesh Parameterization。請參考本站文件「Mesh」。
4. Laplace-Beltrami Operator。請參考本站文件「Triangular Discretization」。
5. Spectral Graph Theory。請參考本站文件「Laplacian Matrix」。

演算法(Jacobi Method)

L是對稱正定矩陣,有唯一最小值,位於x = 0,缺乏討論意義。

兩種解法:一、追加限制條件:已知其中幾點的x值。二、投影:刪去最小的特徵值(零)的特徵向量。

                     vextex label is x                   mesh
min dirichlet energy =================> min sum ‖f′(x)‖² ====>
                   laplacian          differentiate x
min Wij (xi - xj)² ========> min x*Lx ===============> 
                  L is symmetrix
solve (L*+L)x = 0 ==============> solve 2Lx = 0 ===> solve Lx = 0

solve Lx = 0 by gauss-seidel method  <===>  iterative smoothing by average
smallest eigenvectors: min x*Lx             ---> solve Lx = λx
laplacian eigenmaps:   min x*Lx st x*Dx = 1 ---> solve Lx = λDx
best method:           min x*Lx with bc     ---> solve Lx = 0 some xi are known

min x*Lx with bc ---> some xi are known ---> constant b ---> solve Lnew x = b
---> Lnew is block laplacian, rank is low, incosistant system
https://math.stackexchange.com/questions/1727204/

min x*Lx with bc ---> some xi are known ---> solve Lnew x = b by fp iteration
---> solve Lx = 0 by fp iteration ---> iterative smoothing by neighbor
---> what if use fp iteration to solve inconsistent system?

solve Lx = λDx ---> solve D⁻¹Lx = λx
---> D is weight sum and diagonal ---> D⁻¹L is row normalization of L
---> best solution is eigenvector [1,1,1,1,1,1] eigenvalue 0 (useless)
---> 2nd solution as x-coord, 3nd solution as y-coord (good approximation)
---> instead of D⁻¹L, why not L?

Quadratic Form進階變化

二次型除以x長度平方,或者說,二次型的x是單位向量,恰是「虛擬特徵值Rayleigh Quotient」。

二次型等於定值,或者說,二次型的等高線,恰是「圓錐曲線Conic Section」。

二次型大於零,或者說,二次型是正定矩陣,形成「橢圓拋物面函數Elliptic Paraboloid Function」。

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

二次型是對稱正定矩陣,或者說,向量變換後的長度平方,形成「核函數Kernel Function」。

長度函數,乘上正倍率、相加,仍是長度函數。

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 semidefinite)
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))

Quadratic Programming

Quadratic Programming

「二次規劃」。NP-hard。大家轉為討論特例:正定矩陣。

請參考維基百科:

https://en.wikipedia.org/wiki/Quadratic_programming

Semidefinite Programming

Semidefinite Programming

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

Convex Programming

Convex Programming

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

線性規劃、二次規劃、正定規劃,通通是凸規劃的特例。

請參考專著《Convex Optimization》

凸規劃是P問題。然而,判斷凸函數是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