Linear Optimization

Linear Function

「一次函數」。多項式函數,最高次方頂多是一次方。

輸入輸出是數值: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}]

Linear Optimization

「一次最佳化」。一次函數求極值。

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

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

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

‖F(x)‖²有唯一極值,除非遇到特例。

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

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}]]

矩陣微分(矩陣梯度)

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

                              ∂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                             

數學公式(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 Form Optimization

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。兩者的二次型函數圖形一模一樣!因此教科書只討論對稱正定矩陣。

Quadratic Optimization

Quadratic Function

「二次函數」。多項式函數,最高次方頂多是二次方。

輸入輸出是數值: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}]] &)]

二次函數等於定值(二次函數的等高線),恰是「二次曲面Quadric」。兩個變數的情況下,又稱作「圓錐曲線Conic Section」。

g = Plot3D[x*x-y*y+x*y, {x, -3, 3}, {y, -3, 3}, PlotRange -> {-3, 3}, BoxRatios->{1,1,1}, ClippingStyle -> None, Axes->None, Boxed -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &), MeshStyle->{Red,Thickness[0.01]}, MeshFunctions -> (#3 &), Mesh -> {{1}}, RegionFunction -> Function[{x, y, z}, x*x + y*y + z*z <= 9]];

Quadratic Optimization

「二次最佳化」。二次函數求極值。

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是維度。

演算法(Coordinate Descent Method)

「座標下降法」。朝著座標軸方向,正向逆向皆可,一步走到最低處。然後輪到下一個座標軸。

橢圓軸線與座標軸方位相同:最多N步,必定抵達最小值。
橢圓軸線與座標軸方位不相同:無限多步,逐步趨近最小值。
ContourPlot[3*x*x + y*y + x + y + 1, {x, -3, 3}, {y, -3, 3}, Frame -> False, ContourShading -> None, Contours -> 21]

總而言之,應該要往其他方向走。

演算法(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變號,再套用上圖。

演算法(Multigrid Preconditioned Conjugate Gradient Method)

承上,追加兩個外掛,硬是催出效能:Precondition調整矩陣數值、Multigrid調整矩陣大小。

Precondition是以等量乘法公理,讓原本矩陣乘上自訂矩陣,讓原本矩陣變漂亮、容易求解。自訂矩陣稱作Preconditioner。

Multigrid就是Scaling Method的意思。原本矩陣縮放成各種大小(各種解析度),從小到大,依序計算,整合結果。

詳細內容我沒有學會。大家自己看著辦吧。

可以直接使用matlab的pcg()。

延伸閱讀:Quadratic Function特殊用途

二次函數,矩陣是正定矩陣、負定矩陣,形成「橢圓拋物面函數Elliptic Paraboloid Function」。正定矩陣開口向上,負定矩陣開口向下。開口向上的橢圓拋物面函數,乘上正倍率、相加,仍是開口向上的橢圓拋物面函數。

g = Plot3D[(x*x+y*y+x*y-x-y-1)+0*(2*x*x+y*y+2*x*y+x), {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}]] &)];

二階泰勒展開式是二次函數,矩陣是Hessian Matrix。根據Hessian Matrix是否正定,可以初步推測該處是極小值還是鞍點。

syms x y f = y*exp(x - 1) - x*log(y); taylor(f, [x, y], [1, 1], 'Order', 3)
g1 = Plot3D[y*Exp[x-1] - x*Log[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}]] &)]; g2 = Plot3D[x + ((x-1)^2)/2 + ((y-1)^2)/2, {x, -3, 3}, {y, -3, 3}, PlotRange -> {-3, 3}, BoxRatios->{1,1,1}, ClippingStyle -> None, Mesh-> None, Axes->None, Boxed -> False, ColorFunction -> "SolarColors"];

延伸閱讀:Quadratic Form特殊用途

二次型除以x長度平方,換句話說,向量經過線性變換再投影回去,即是「虛擬特徵值Rayleigh Quotient」。

二次型開根號(必須是半正定矩陣),換句話說,向量經過線性變換的長度,形成「歐氏長度函數Euclidean Length Function」。

Linear Programming(Under Construction!)

Programming

「規劃」是指約束最佳化。已經發展出許多類型:

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
LP ⊂ QP ⊂ SOCP ⊂ SDP

約束條件,等式通通換成不等式:以大於等於、小於等於的兩道不等式,取代一道等式。

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。大家轉為討論特例:正定矩陣。

請參考專著《Numerical Optimization》

Integer Programming

Integer Programming / Mixed Integer Programming

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

「混和整數規劃」。有些變數不必是整數。

工具如GurobiCPLEX

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

Branch and Bound(分支定界法)

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

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

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

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

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

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

Convex Optimization(Under Construction!)

Convex Optimization

「凸最佳化」。目標函數是凸函數。

一次最佳化、二次最佳化,通通是凸最佳化的特例。

請參考專著《Convex Optimization》

Convergence Rate

http://www.yisongyue.com/courses/cs155/2018_winter/lectures/Lecture_02.pdf

f is convex and f(x) - f(x*) ≤ ε
|f(a) - f(b)| ≤ ρ |a - b|                    O(1/ε²)  [Lipschitz]
|∇f(a) - ∇f(b)| ≤ ρ |a - b|                 O(1/ε)
f(a) ≥ f(b) + ∇f(b)ᵀ(a-b) + 1/2 ρ (a-b)²    O(log 1/ε)  [Taylor 2nd]

Convex Programming(Under Construction!)

Convex Programming

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

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

凸規劃是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

Duality in Optimization(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 cᵀx   subject to Ax ≤ b
min bᵀy   subject to Aᵀy ≥ c
Ax ≤ b  ---> yᵀAx ≤ yᵀb     (y ≥ 0)
Aᵀy ≥ c ---> xᵀAᵀy ≥ xᵀc    (x ≥ 0)

cᵀx = xᵀc ≤ xᵀAᵀy = yᵀAx ≤ yᵀb = bᵀy
Complementary Slackness Theorem
https://www.linearprogramming.info/primal-dual-relationships-in-linear-programming-duality-theory-in-lp/
https://s3.amazonaws.com/content.udacity-data.com/courses/gt-cs6505/duality.html

max 1 x1 + 2 x2              (x1 >= 0, x2 >= 0)
{ 3 x1 + 4 x2 <= 7
{ 5 x1 + 6 x2 <= 8

{ (3 x1 + 4 x2) y1 <= 7 y1   (y1 >= 0, y2 >= 0)
{ (5 x1 + 6 x2) y2 <= 8 y2

3 x1y1 + 4 x2y1 + 5 x1y2 + 6 x2y2 <= 7 y1 + 8 y2

(3 y1 + 5 y2) x1 + (4 y1 + 6 y2) x2 <= 7 y1 + 8 y2s

min 7 y1 + 8 y2
{ (3 y1 + 5 y2) >= 1
{ (4 y1 + 6 y2) >= 2