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

二次函數,矩陣是正定矩陣、負定矩陣,稱作「橢圓拋物面函數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}]] &)];

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。

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

當A是正定矩陣、負定矩陣,「一次微分等於零」的地方必是極值。f′(x) = (A + Aᵀ)x + b = 0,解x,就是極值。

問題化作矩陣求解(A + Aᵀ)x = -b。高斯消去法O(N³)。鬆弛法O(N²T)。N是維度。T是遞推次數。

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

延伸閱讀:數學證明

{ 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(dᵢ, Adⱼ) = 0:最多N步,必定抵達最小值。

「共軛方向法」是第三種方式。

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

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

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

延伸閱讀:數學證明

一、最小值位置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

四、任意向量v的性質:前進方向向量互相線性獨立,N項必能湊得任意向量。
v = c₀d₀ + c₁d₁ + ... + cN-1dN-1

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

六、步伐大小a的性質:前進方向互相傍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₂, ...

Conjugate Gradient Method(共軛梯度法)

承上,以遞推公式產生每一步的前進方向。

本次梯度反方向(朝向低處)、上次前進方向,兩者的加權平均值,做為本次前進方向。

亦可解釋成:N個前進方向,是將每一處的梯度反方向-f′(x)扳成互相傍A垂直,邊走邊扳。

共軛梯度法的時間複雜度也是O(N³),但是計算步驟更少。

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

上圖當中,A是對稱正定矩陣,f(x) = 1/2 xᵀAx - bᵀx,f′(x) = Ax - b。r = -f′(x)。p = d。α = step。

如果要最佳化f(x) = xᵀAx + bᵀx + c,記得先將矩陣化做對稱正定矩陣A、將A乘以2、將b變號,再套用上圖。

延伸閱讀:數學證明

==================
quadratic function
==================

A is symmetric positive definite
f(x) = 1/2 xᵀAx - bᵀx     // people use this formulation
f′(x) = Ax - b            // because its gradient is beautiful

=========
algorithm
=========

x0 = random vector                 // initial position
g0 = f′(x0) = A x0 - b             // initial gradient
d0 = -g0                           // initial direction
if (|g0| < ε) return x0
for (k = 0; k < n; k++)            // loop all dimensions
{
    αk = ❮gk, gk❯ / ❮dk, A dk❯     // update step size
    xk+1 = xk + αk dk              // update position
    gk+1 = gk + αk A dk            // update gradient
    if (|gk| < ε) return xk+1
    βk = ❮gk+1, gk+1❯ / ❮gk, gk❯   // update direction coefficient
    dk+1 = -gk+1 + βk dk           // update direction
}
return xk+1

================
steepest descent
================

1. next position: xk+1
xk+1 = xk + αk dk       // step from xk to xk+1 among direction dk

2. next gradient: gk+1
gk+1 = f′(xk+1)
gk+1 = f′(xk + αk dk)
gk+1 = A (xk + αk dk) - b
gk+1 = (A xk - b) + αk A dk
gk+1 = gk + αk A dk        // iteration formula

3. this direction and next gradient are orthogonal
min f(xk+1)
min f(xk + αk dk)
df/dαk = 0
f′(xk + αk dk)ᵀ dk = 0
f′(xk+1)ᵀ dk = 0
gk+1ᵀ dk = 0               // dk is the tangent line to levelset
❮gk+1, dk❯ = 0             // xk+1 is the tangent point

4. step size: αk   [incorrect result]
f′(xk + αk dk)ᵀ dk = 0
[f′(xk) + s f′(dk)]ᵀ dk = 0          // f′(x) is linear(calculus)
αk = - ❮gk, dk❯ / ❮f′(dk), dk❯       // but not linear(algebra)

4. step size: αk
{ gk+1ᵀ dk = 0             // tangent
{ gk+1 = gk + αk A dk      // iteration formula
dkᵀ gk+1 = dkᵀ gk + αk dkᵀ A dk
       0 = dkᵀ gk + αk dkᵀ A dk
αk = - ❮dk, gk❯ / ❮dk, A dk❯

5. gradient delta = A(position delta)   [no use]
xk+1 = xk + αk dk        |   gk+1 = gk + αk A dk
xk+1 - xk = αk dk        |   gk+1 = gk + A(xk+1 - xk)
A(xk+1 - xk) = αk A dk   |   gk+1 - gk = A(xk+1 - xk)

6. A(optimal solution) = b
min f(x)
f′(x*) = 0
Ax* - b = 0
Ax* = b

7. gradient = A(error)
gk = f′(xk)
gk = A xk - b
gk = A xk - A x* 
gk = A (xk - x*)         // people think of A-orthogonal

==========================
conjugate direction method
==========================

1. A-othrogonal directions: dk
❮di, A dj❯ = 0

2. dot product provides coefficient
assume x* = c0 d0 + c1 d1 + ...
diᵀ A x* = ci diᵀ A di
diᵀ b    = ci diᵀ A di
ci = diᵀ b / diᵀ A di
ci = ❮di, b❯ / ❮di, A di❯

3. A-othrogonal directions are linear independent
proof by contradiction
assume dk = c0 d0 + c1 d1 + ... + ck-1 dk-1 is not linear independent
{ dkᵀ A dk > 0     // A is symmetric positive definite and dk != 0
{ dkᵀ A dk = c0 dkᵀ A d0 + c1 dkᵀ A d1 + ... + ck-1 dkᵀ A dk-1 = 0

4. n step ensures optimal solution
x* - x0 = α1 d1 + α2 d2 + ... + αn-1 dn-1

5. position iteration
xk+1 = xk + αk dk
x1 - x0 = α0 d0
x2 - x0 = α0 d0 + α1 d1
:    :      :       :
x* - x0 = α0 d0 + α1 d1 + ... + αn-1 dn-1

6. same step size! conjugate direction method is steepest descent!
αk = dkᵀ A (x* - x0) / dkᵀ A dk
αk = dkᵀ A (x* - xk) / dkᵀ A dk
αk = - dkᵀ gk / dkᵀ A dk
αk = - ❮dk, gk❯ / ❮dk, A dk❯

7. gradient iteration
gk+1 = gk + αk A dk
g1 = g0 + α0 A d0
g2 = g0 + α0 A d0 + α1 A d1
:    :       :         :
gk = g0 + α0 A d0 + α1 A d1 + ... + αk-1 A dk-1

8. alternative gradient
{ dkᵀ g0~k     = dkᵀ g0
{ dkᵀ gk+1~n-1 = dkᵀ g0 + αk dkᵀ A dk = 0    // dkᵀ gk+1 = 0
{ diᵀ gj = diᵀ g0  when i >= j
{ diᵀ gj = 0       when i < j                // dkᵀ gk+1 = 0

9. alternative step size: αk   [no use]
dkᵀ gk+1 = 0
dkᵀ gk+1 = dk ᵀ g0 + αk dkᵀ A dk = 0
αk = - ❮dk, g0❯ / ❮dk, A dk❯

10. Dk = span(d1, d2, ..., dk) and next gradient are orthrogonal
diᵀ gk+1 = 0   when i < k+1

11. first k step provides optimal solution of "min f(x) when x in Dk"
min f(x0 + sum ck dk)        |   let ck = αk
df/dck = 0                   |   f′(x0 + sum ck dk) = gk+1
f′(x0 + sum ck dk)ᵀ dk = 0   |   gk+1ᵀ dk = 0

=========================
conjugate gradient method
=========================

1. next direction: dk+1
dk+1 = -gk+1 + βk dk           // negative gradient goto minimum

2. direction coefficient: βk   [however dk+1 is unknown]
{ dk+1 = -gk+1 + βk dk         // adaptive gradient descent
{ dkᵀ gk+1 = 0                 // steepest descent
dk+1 + gk+1 = βk dk
dkᵀ dk+1 + dkᵀ gk+1 = βk dkᵀ dk
dkᵀ dk+1 = βk dkᵀ dk
βk = ❮dk, dk+1❯ / ❮dk, dk❯

3. direction coefficient: βk   [however precision is not good?]
{ dk+1 = -gk+1 + βk dk         // adaptive gradient descent
{ dk+1ᵀ A dk = 0               // conjugate direction method
(-gk+1 + t dk)ᵀ A dk = 0
-gk+1ᵀ A dk + βk dkᵀ A dk = 0
βk = ❮gk+1, A dk❯ / ❮dk, A dk❯

4. all gradients are orthogonal
diᵀ gj = 0   when i < j
{ (-gi + ti-1 di-1)ᵀ gj = 0   when i > 0 and i < j
{ let di = -gi                when i = 0 and i < j
giᵀ gj = 0   when i < j

5. initial direction: d0 = -g0
such that all gradients are orthogonal

6. direction -> gradient
{ dk+1 = -gk+1 + βk dk         // adaptive gradient descent
{ dkᵀ gk+1 = 0                 // steepest descent
dk+1ᵀ gk+1 = -gk+1ᵀ gk+1

7. alternative gradient
{ dkᵀ g0~k     = dkᵀ g0
{ dkᵀ gk+1~n-1 = dkᵀ g0 + αk dkᵀ A dk = 0
dk+1ᵀ g0~k+1 = -gk+1ᵀ gk+1

8. alternative direction coefficient: βk
dk+1 = -gk+1 + βk dk
dk+1ᵀ gk = -gk+1ᵀ gk + βk dkᵀ gk  // 1. alternative gradient
-gk+1ᵀ gk+1 = 0 - βk gkᵀ gk       // 2. all gradients are orthogonal
βk = ❮gk+1, gk+1❯ / ❮gk, gk❯      // 3. direction -> gradient

9. alternative step size: αk
αk = - ❮dk, gk❯ / ❮dk, A dk❯
αk = ❮gk, gk❯ / ❮dk, A dk❯

Multigrid Preconditioned Conjugate Gradient Method(多重網格預條件共軛梯度法)

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

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

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

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

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

Convex Optimization

Convex Optimization

「凸最佳化」。凸函數求極值。

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

演算法設計:大家沿用一般函數的演算法,例如梯度下降法。雖然凸函數有專門演算法,但是不實用,例如橢球法。

演算法分析:一般函數,實施演算法,不保證趨近極值位置(收斂),無法推導時間複雜度。但是凸函數可以。

Rate of Convergence

https://www.cs.ubc.ca/~schmidtm/Courses/540-W19/L5.pdf

誤差:當前位置與極值位置的差距。

收斂速率:走一步,前後兩處的誤差比率。

各種演算法、各種特殊函數,收斂速率皆不同。利用收斂速率,判斷時間複雜度:達到預定的數值精確度,需要走多少步。

例如梯度下降法用於凸函數,時間複雜度為O(log 1/ε)。

run algorithm until |f(x) - f(x*)| ≤ ε

f is Lipschitz:  |f(a) - f(b)| ≤ ρ |a - b|
O(1/ε²) iterations

∇f is Lipschitz: |∇f(a) - ∇f(b)| ≤ ρ |a - b|
O(1/ε) iterations

f is strong convex: f(a) ≥ f(b) + ∇f(b)ᵀ(a-b) + 1/2 ρ (a-b)²
O(log 1/ε) iterations

Disciplined Convex Optimization

凸最佳化是P問題。然而,判斷凸函數是NP-complete問題。

於是有人改弦易轍,排列組合常見運算符號、常見函數,藉由已知的數學性質,以判斷凸函數,甚至判斷要使用哪個演算法。

請參考專論《Disciplined Convex Programming》

照理來說,這是在判斷凸函數,不應該取名為最佳化。

Ellipsoid Method(橢球法)

https://stanford.edu/class/ee364b/lectures/ellipsoid_method_notes.pdf

製造一個橢球,涵蓋極值位置。依照當前梯度,移動橢球球心、調整橢球形狀、縮小橢球體積,球心逐步收斂至極值位置。

僅有理論上的價值,沒有實務上的價值。

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

現代社會經常使用數學規劃,例如經濟交易、交通運輸、工業生產。數學規劃儘管不是電資科系的學習重點,然而卻是商管科系的必修專業。

經典數學領域亦可使用數學規劃,例如組合最佳化、排程理論。數學規劃比起經典組合算法,有過之而無不及,計算時間更短。

知名工具如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)個頂點。

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)

https://ocw.mit.edu/courses/sloan-school-of-management/15-053-optimization-methods-in-management-science-spring-2013/tutorials/MIT15_053S13_tut07.pdf

實施消去法,極少數情況下,最大值持續不變,不斷鬼打牆。

計算順序(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

Quadratic Programming(Under Construction!)

Quadratic Programming

「二次規劃」。目標函數是二次函數,約束函數是一次函數。

NP-hard。大家轉為討論正定矩陣。

大家沿用內點法。

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,求解。兩解相減做為行進方向,步伐大小盡量大,更換限制條件。不斷行進,直到乘數非負。

二次函數,偏微分等於零,即是一次方程組,可以用高斯消去法求解。

Wolfe's Modified Simplex Method(單形法)

二次規劃化作一次規劃。

拉格朗日乘數、單形法、限制某些變數不可同時非零。

請見專著《Operations Research: Applications and Algorithms》。

Sequential Quadratic Programming(循序二次規劃)
Lagrange-Newton Method(拉格朗日牛頓法)

約束最佳化化作二次規劃。

請見專著《Numerical Optimization》

Trust Region Method(信賴區域法)

最佳化化作二次最佳化。

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

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

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

Convex Programming

Convex Programming

「凸規劃」。目標函數是凸函數,可行解是凸集合。

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

請見專著《Convex Optimization》

Integer Programming

Integer Programming / Mixed Integer Programming

「整數規劃」。所有變數必須是整數。整數解。

「混合整數規劃」。有些變數不必是整數。部分整數解。

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

請見專著《Operations Research: Applications and Algorithms》。

Branch and Bound(分支定界法)

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

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

Cutting Plane Method(切割平面法)

整數一次規劃:分數解所在的約束條件,任取其中一個,擷取小數部分,使之小於1(小於等於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

Fenchel 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://mdav.ece.gatech.edu/ece-8823-spring2019/notes/06-fenchel-dual.pdf
https://people.ok.ubc.ca/bauschke/Research/68.pdf
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(Complementary Slackness)

一次規劃的對偶。

原始問題稱作primal,對偶問題稱作dual,兩個問題等價。根據單形法作者的回憶錄,primal是他爹創造的名稱。

primal: max cᵀx   subject to Ax ≤ b, x ≥ 0
dual:   min bᵀy   subject to Aᵀy ≥ c, y ≥ 0
primal:                       dual:
-------------------------     |   ∑     ∑     ∑   |  ∑
∑ a₁₁x₁ a₁₂x₂ a₁₃x₃ ≤ b₁      | a₁₁y₁ a₁₂y₁ a₁₃y₁ | b₁y₁
∑ a₂₁x₁ a₂₂x₂ a₂₃x₃ ≤ b₂      | a₂₁y₂ a₂₂y₂ a₂₃y₂ | b₂y₂
∑ a₃₁x₁ a₃₂x₂ a₃₃x₃ ≤ b₃      | a₃₁y₃ a₃₂y₃ a₃₃y₃ | b₃y₃
∑ a₄₁x₁ a₄₂x₂ a₄₃x₃ ≤ b₄      | a₄₁y₄ a₄₂y₄ a₄₃y₄ | b₄y₄
-------------------------     |   |V    |V    |V  |  ||
∑  c₁x₁  c₂x₂  c₃x₃ = max     |   c₁    c₂    c₃  |  min

一次規劃的對偶,源自窮舉法。

一、以約束條件夾擠最大值,得到對偶問題的目標函數。

約束條件的加權總和,不影響正解。窮舉權重,約束條件的加權總和與目標函數,取其一致者;常數項,取其最小者。

二、最大值不超過約束條件,得到對偶問題的不等式。

約束條件的加權總和,重新整理成變數們的加權總和。窮舉權重,並讓第一個變數以外的權重是零,得到第一道不等式。

"min bᵀy"
row₁ k₁ + row₂ k₂ + ...       ≤ b₁k₁ + b₂k₂ + ...   row sum, kᵢ ≥ 0
∑ kᵢaᵢ₁ x₁ + ∑ kᵢaᵢ₂ x₂ + ... ≤ b₁k₁ + b₂k₂ + ...   column sum
c₁x₁ + c₂x₂ + ...             ≤ b₁k₁ + b₂k₂ + ...   ∑ kᵢaᵢ₁ = c₁ for some k
objective                     = min bᵀy             k is y
"subject to Aᵀy ≥ c"
cᵀx               ≤ row₁ k₁ + row₂ k₂ + ...         row sum
c₁x₁ + c₂x₂ + ... ≤ ∑ kᵢaᵢ₁ x₁ + ∑ kᵢaᵢ₂ x₂ + ...   column sum
c₁x₁              ≤ ∑ kᵢaᵢ₁ x₁                      ∑ kᵢaᵢ₂ = c₂ for some k
c₁                ≤ ∑ kᵢaᵢ₁                         divide x₁ when x₁ ≠ 0
c                 ≤ Aᵀy                             k is y

一次規劃的對偶,本質是拉格朗日對偶。

max cᵀx   subject to Ax ≤ b, x ≥ 0
L(x,y) = cᵀx - yᵀ(Ax - b)
L(x,y) = cᵀx - yᵀAx + yᵀb
L(x,y) = bᵀy - xᵀAᵀy + xᵀc
L(x,y) = bᵀy - xᵀ(Aᵀy - c)
min bᵀy   subject to Aᵀy ≥ c, y ≥ 0

原始問題的最佳解等於對偶問題的最佳解。

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

一次規劃是特例,一次不等式是通例。

原始問題的約束條件、對偶問題的約束條件、兩個問題的目標函數相等,三者聯立得到一次不等式。

linear programming          linear inequality
maximize   cᵀx              { Ax ≤ b    (primal constraint)
subject to Ax ≤ b    --->   { Aᵀy ≥ c   (dual constraint)
           x ≥ 0            { cᵀx = bᵀy (two objectives equality)