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}]]
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                    移項,當A的rank足夠時

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是遞推次數。

正定矩陣A化作對稱正定矩陣(A + Aᵀ) / 2,二次函數圖形一模一樣,極值一模一樣。因此大家習慣直接討論對稱正定矩陣。

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

問題化作對稱正定矩陣求解(2A)x = -b。Cholesky分解O(N³)。鬆弛法O(N²T)。N是維度。T是遞推次數。

linear equation: conjugate gradient method

人為加工

二次最佳化的主要用途是一次方程組求解。當A是對稱正定矩陣,f(x) = 1/2 xᵀAx - bᵀx的最小值位置,即是Ax = b的解。

(正定確保唯一解,對稱確保形成A。)

             a unique solution exist
argmin f(x) <=======================> solve f′(x) = 0

let f(x)  = 1/2 xᵀAx - bᵀx
    f′(x) = 1/2 (A+Aᵀ)x - b

     argmin 1/2 xᵀAx - bᵀx
---> solve 1/2 (A+Aᵀ)x - b = 0   (A is positive definite)
---> solve Ax = b                (A is symmetric)

這導致世上所有的二次最佳化演算法,全部都是針對f(x) = 1/2 xᵀAx - bᵀx,而非f(x) = xᵀAx + bᵀx + c。非常奇葩。

想要最佳化標準的二次函數f(x) = xᵀAx + bᵀx + c,記得先將A乘以2、將b變號、將c刪掉,才能實施下述演算法。

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]

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

line search(直線搜尋)

給定地點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]

也就是說,最低處的梯度方向f′(xnext) = Axnext - b、行進方向d,互相垂直,點積是零。

藉此推導步伐大小α = ❮-f′(x) , d❯ / ❮Ad , d❯。

數學公式:

α = ❮-f′(x) , d❯ / ❮Ad , d❯
  = ❮b - Ax , d❯ / ❮Ad , d❯
xnext = x + α d
d: direction
α: step size
f′(x): gradient at position x        (normally ∇f(x))
❮a,b❯: dot product of vector a and b (normally a∙b)

演算法步驟:

x = random vector
while (true)
{
    d = get_direction(x);
    α = ❮b - Ax , d❯ / ❮Ad , d❯
    if (α < ε) return x
    x = x + α d
}

延伸閱讀:推導過程

f(x)  = 1/2 xᵀAx - bᵀx
f′(x) = Ax - b
⎰ xnext = x + α d
⎱ α = argmin f(x + α d)

❮f′(xnext)      , d❯ = 0          perpendicularity
❮f′(x + α d)    , d❯ = 0          substitution
❮A(x + α d) - b , d❯ = 0          substitution
❮Ax + α Ad  - b , d❯ = 0          expansion
❮Ax - b + α Ad  , d❯ = 0          arrangement
❮Ax - b , d❯ + α ❮Ad , d❯ = 0     distributive law
α = - ❮Ax - b , d❯ / ❮Ad , d❯     transposition
α =   ❮b - Ax , d❯ / ❮Ad , d❯     arrangement
α =   ❮-f′(x) , d❯ / ❮Ad , d❯     just a coincidence

steepest descent method(最陡下降法)

承上,仔細規定每一步的前進方向:梯度反方向(朝向低處)。

d = -f′(x) = b - Ax
α = ❮d , d❯ / ❮Ad , d❯
xnext = x + α d

由於人為加工,梯度反方向-f′(x) = b - Ax碰巧等於殘差r = b - Ax。有些書籍將前進方向d改寫成殘差r,看上去比較帥氣。

r = b - Ax
α = ❮r , r❯ / ❮Ar , r❯
xnext = x + α r
r: residual (equal to negative gradient. just a coincidence.)

conjugate direction method(共軛方向法)

承上,仔細規定每一步的前進方向:互相傍A垂直。

一、前後垂直  ❮dᵢ , dᵢ₊₁❯ = 0:逐步趨近最小值。
二、互相垂直  ❮dᵢ , dⱼ  ❯ = 0:最多N步,通常無法抵達最小值。
三、互相傍A垂直❮dᵢ , Adⱼ ❯ = 0:最多N步,必定抵達最小值。

第一種方式是「最陡下降法」。第三種方式是「共軛方向法」。

挑選N個前進方向,互相傍A垂直,有無限多種方式。步伐大小符合直線搜尋,就能保證N步抵達最小值。

挑選N個前進方向,一種方式是Krylov subspace聯合Gram–Schmidt orthonormalization,將v, Av, A²v, A³v, ..., Aᴺ⁻¹v扳成互相傍A垂直,邊走邊扳。v可以是任意向量,實務上是b。

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

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

數學公式:

α = ❮d , -g❯ / ❮d , A d❯
xnext = x + α d
gnext = g + α A d
dnext = A-orthonomalize(A d)
g: gradient  f′(x) = Ax - b
r: residual  -f′(x) = b - Ax

演算法步驟:

x₀ = random vector                // initial position
g₀ = f′(x₀) = Ax₀ - b             // initial gradient
d₀ = random vector                // initial direction
d₀ = d₀ / ‖d₀‖                    // normalization
if (‖g₀‖ < ε) return x₀           // ‖g₀‖ = sqrt(❮g₀, g₀❯)
for (k = 0 ... n-1)               // loop all dimensions
{
    αk = ❮dk, -gk❯ / ❮dk, A dk❯   // update step size
    xk+1 = xk + αk dk             // update position
    gk+1 = gk + αk A dk           // update gradient
    if (‖gk+1‖ < ε) return xk+1   // ‖gk+1‖ = sqrt(❮gk+1, gk+1❯)
    dk+1 = A dk                   // update direction
    for (t = 0 ... k)             // A-orthonormalization
        dk+1 -= ❮dk+1, A dt❯ dt  
    dk+1 = dk+1 / ‖dk+1‖          // normalization
}
return xk+1

延伸閱讀:推導過程

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

二、誤差e = x* - x的性質:誤差經過A變換,碰巧是梯度反方向。
Ae = Ax* - Ax = (Ax* - b) - (Ax - b) = f′(x*) - f′(x) = -f′(x)

三、前進方向d的性質:前進方向互相傍A垂直。
❮dᵢ, Adⱼ❯ = 0

四、前進方向d的性質:若前進方向互相傍A垂直,則線性獨立。
A-orthogonal directions are linear independent
proof by contradiction
assume dₖ = c₀d₀ + c₁d₁ + ... + cₖ₋₁dₖ₋₁ is not linear independent
1. ❮dₖ , Adₖ❯ = c₀ ❮d₀ , Adₖ❯ + ... + cₖ₋₁ ❮dₖ₋₁ , Adₖ❯ = 0
   (A-orthogonal directions)
2. ❮dₖ , Adₖ❯ > 0
   (A is symmetric positive definite and dₖ ≠ 0)

五、任意向量v的性質:若前進方向互相線性獨立,則N步必能形成任意向量。
v = c₀d₀ + c₁d₁ + ... + cɴ₋₁dɴ₋₁

六、初始誤差e₀ = x* - x₀的性質:N步必能走到最佳解。步伐大小α。
e₀ = α₀d₀ + α₁d₁ + ... + αɴ₋₁dɴ₋₁

七、步伐大小α的性質:恰好等同直線搜尋。
e₀ = α₀d₀ + α₁d₁ + ... + αɴ₋₁dɴ₋₁
Ae₀ = A(α₀d₀ + α₁d₁ + ...)
❮d₀ , Ae₀❯ = ❮d₀ , A(α₀d₀ + α₁d₁ + ...)❯
❮d₀ , Ae₀❯ = α₀ ❮d₀ , Ad₀❯ + α₁ ❮d₀, Ad₁❯ + ...
❮d₀ , Ae₀❯ = α₀ ❮d₀ , Ad₀❯
α₀ = ❮d₀, Ae₀❯ / ❮d₀, Ad₀❯
α₀ = ❮d₀, -f′(x₀)❯ / ❮d₀, Ad₀❯     line search

similarly,
consider Ae₁ and get α₁ = ❮d₁, -f′(x₁)❯ / ❮d₁, Ad₁❯
consider Ae₂ and get α₂ = ❮d₂, -f′(x₂)❯ / ❮d₂, Ad₂❯
   :

八、N處梯度g的性質:梯度遞推公式gnext = g + α A d。
g₁ - g₀ = f′(x₁) - f′(x₀) = -Ae₁ + Ae₀ = A(-e₁ + e₀)
        = A(x₁ - x* + x* - x₀) = A(x₁ - x₀) = A(α₀ d₀) = α₀ A d₀
g₁ = g₀ + α₀ A d₀
g₂ = g₀ + α₀ A d₀ + α₁ A d₁
g₃ = g₀ + α₀ A d₀ + α₁ A d₁ + α₂ A d₂
   :

九、N處位置x的性質:第k步走到當前最佳解(考慮第k個方向能夠觸及的所有位置)。
  N處位置x的性質:頭k步走到當前最佳解(考慮頭k個方向能夠觸及的所有位置)。
  證明省略。
xk is the optimal solution of min f(x) when x ∈ Dk
D₀ = span(d₀) = span(v)
x₁ = x₀ + α₀d₀ = argmin { f(x) : x ∈ D₀ }
D₁ = span(d₀, d₁) = span(v, Av)
x₂ = x₀ + α₀d₀ + α₁d₁ = argmin { f(x) : x ∈ D₁ }
D₂ = span(d₀, d₁, d₂) = span(v, Av, A²v)
x₃ = x₀ + α₀d₀ + α₁d₁ + α₂d₂ = argmin { f(x) : x ∈ D₂ }
   :

十、N處梯度g的性質:垂直於先前所有方向。
  證明省略。
g₁ ⟂ D₀ = span(d₀)
g₂ ⟂ D₁ = span(d₀, d₁)
g₃ ⟂ D₂ = span(d₀, d₁, d₂)
        :

conjugate gradient method(共軛梯度法)

承上,仔細規定拿來扳正的向量:梯度反方向。

一、任意向量v。v, Av, A²v, A³v, ..., Aᴺ⁻¹v,扳成互相傍A垂直。
二、N個位置的梯度反方向-f′(x₀) ... -f(xɴ₋₁),扳成互相傍A垂直。

第一種方式是「共軛方向法」。第二種方式是「共軛梯度法」。

邊走邊扳。利用遞推公式,製造每一步的前進方向:本次梯度反方向、上次前進方向,兩者的加權平均數,作為本次前進方向。

邊走邊扳。N個位置的梯度反方向,扳正之後所得到的N個行進方向,恰巧導致這N個梯度互相垂直!遞推公式可以稍作簡化!

也就是說,共軛梯度法同時滿足兩個性質:

一、前進方向互相傍A垂直❮dᵢ , Adⱼ❯ = 0
二、梯度互相垂直❮gᵢ , gⱼ❯ = 0
 (殘差互相垂直❮rᵢ , rⱼ❯ = 0)

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

數學公式:

α = ❮g, g❯ / ❮d, A d❯
xnext = x + α d
gnext = g + α A d
β = ❮gnext, gnext❯ / ❮g, g❯
dnext = -gnext + β d

演算法步驟:

x₀ = random vector                 // initial position
g₀ = f′(x₀) = Ax₀ - b              // initial gradient
d₀ = -g₀                           // initial direction
if (‖g₀‖ < ε) return x₀            // ‖g₀‖ = sqrt(❮g₀, g₀❯)
for (k = 0 ... n-1)                // 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+1‖ < ε) return xk+1    // ‖gk+1‖ = sqrt(❮gk+1, gk+1❯)
    βk = ❮gk+1, gk+1❯ / ❮gk, gk❯   // update direction coefficient
    dk+1 = -gk+1 + βk dk           // update direction
}
return xk+1

演算法步驟(殘差風格),借用一下維基百科的圖片:

延伸閱讀:推導過程

1. next direction: dk+1
⎰ d0   = -g0                    // negative gradient goto minimum
⎱ dk+1 = -gk+1 + βk dk          // adaptive gradient descent

2. a special formula: ❮dk , gk❯ = ❮-gk , gk❯
⎰ dk+1 = -gk+1 + βk dk          // adaptive gradient descent
⎱ ❮dk , gk+1❯ = 0               // line search
❮dk+1 , gk+1❯ = ❮-gk+1 , gk+1❯ + βk ❮dk , gk+1❯
❮dk+1 , gk+1❯ = ❮-gk+1 , gk+1❯
❮dk , gk❯ = ❮-gk , gk❯          // dk = -gk when k = 0

3. alternative step size: αk
αk = ❮dk , -gk❯ / ❮dk , A dk❯   // line search
αk = ❮gk , gk❯  / ❮dk , A dk❯   // 2.

4. direction coefficient: βk
⎰ dk+1 = -gk+1 + βk dk          // adaptive gradient descent
⎱ ❮dk+1 , A dk❯ = 0             // conjugate direction method
❮dk+1 , A dk❯ = ❮-gk+1 , A dk❯ + βk ❮dk , A dk❯
0             = ❮-gk+1 , A dk❯ + βk ❮dk , A dk❯
βk = ❮gk+1 , A dk❯ / ❮dk , A dk❯

5. all gradients are orthogonal
5a. gradient gk+1 is orthogonal to direction d₀ ... dk
❮gk+1 , di❯ = 0   for i = 0 ... k   // conjugate direction method
gk+1 ⟂ Dk = span(d₀, d₁, ..., dk)   // conjugate direction method
5b. take all negative gradients for A-orthonormalization
-gi ∈ Dk = span(d₀, d₁, ..., dk)   for i = 0 ... k
❮gi , gk+1❯ = 0   for i = 0 ... k

6. alternative direction coefficient: βk
βk = ❮gk+1 , A dk❯      / ❮dk , A dk❯
βk = ❮gk+1 , αk A dk❯   / ❮dk , αk A dk❯     // expand a fraction
βk = ❮gk+1 , gk+1 - gk❯ / ❮dk , gk+1 - gk❯   // gradient formula
βk = ❮gk+1 , gk+1❯      / ❮dk , gk+1 - gk❯   // orthogonal gradients
βk = ❮gk+1 , gk+1❯      / ❮dk , -gk❯         // line search
βk = ❮gk+1 , gk+1❯      / ❮gk , gk❯          // 2.

linear equation: biconjugate gradient method🚧

conjugate gradient method(共軛梯度法)

現實世界當中,二次最佳化的主要用途其實是一次方程組求解。首項係數1/2就是為了湊出一次方程組。

     min 1/2 xᵀAx - bᵀx            quadratic optimization
---> solve 1/2 (A+Aᵀ)x - b = 0     極值位於一次微分等於零的地方
---> solve 1/2 (A+Aᵀ)x = b         linear equation
---> solve Ax = b                  if A = Aᵀ

當矩陣是正定矩陣,導致二次函數存在唯一極值,導致一次函數存在唯一解。兩者答案相同。

當矩陣是對稱矩陣(故為對稱正定矩陣),一次函數的矩陣是A。實施二次最佳化演算法(例如共軛梯度法),就能求解。

對於對稱正定矩陣,數學家採用共軛梯度法,捨棄Cholesky分解。

biconjugate gradient method(雙共軛梯度法)

     min 1/2 xᵀAx - bᵀx            quadratic optimization
---> solve 1/2 (A+Aᵀ)x - b = 0     極值位於一次微分等於零的地方
---> solve 1/2 (A+Aᵀ)x = b         linear equation

當矩陣不是對稱矩陣(只是一般的正定矩陣),一次方程組的矩陣是(A+Aᵀ)/2而不是A。

(A+Aᵀ)/2是對稱矩陣。想盡辦法調整A,一次方程組的矩陣仍是對稱矩陣,無法得到非對稱矩陣。

對於非對稱的正定矩陣,數學家硬是將共軛梯度法魔改成雙共軛梯度法,硬是將(A+Aᵀ)/2魔改成A。然而,演算法失去了原本的數學性質,不見得逐步趨近正解。演算法可能永不結束。求解必須拚運氣,看準時機脫手離場。

由於不見得趨近正解,也沒必要執著於正定矩陣了。乾脆處理任意矩陣吧。

stabilized biconjugate gradient method(穩定雙共軛梯度法)

conjugate gradient squared method(共軛梯度平方法)

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

linear equation: “genuine” minimal residual method

minimal residual method(最小殘差法)

《Annotated Algorithms in Python: With applications in Physics, Biology, and Finance》

原版的最小殘差法,用途更加刁鑽、過程更加複雜──求得對稱矩陣的最小平方誤差解。

此處的最小殘差法,則是該書作者重新精煉而得,名稱與內容更加相符。我也認為這才是名副其實的最小殘差法。

最小殘差法用來解Ax = b。已知A和b,求x。

方便起見,假設存在唯一解。

根據反矩陣,Ax = b等價於x = A⁻¹b。

根據power decomposition【尚無正式名稱】,A⁻¹是A⁰ ... Aᴺ⁻¹的線性組合(加權總和)。更進一步,x是A⁰b ... Aᴺ⁻¹b的線性組合。更進一步,x是Krylov subspace的線性組合。

想要求x,一種方式是將初始誤差e₀ = x* - x₀投影到Krylov subspace。回憶一下共軛方向法提到的性質:

五、任意向量v的性質:若前進方向互相線性獨立,則N步必能形成任意向量。
v = c₀d₀ + c₁d₁ + ... + cɴ₋₁dɴ₋₁

六、初始誤差e₀ = x* - x₀的性質:N步必能走到最佳解。步伐大小α。
e₀ = α₀d₀ + α₁d₁ + ... + αɴ₋₁dɴ₋₁

理想方式是平行投影,而平行投影就是求解,事情繞回原點,事情沒有解決。我們只好採用垂直投影,不保證N步抵達最小值,也不保證逐步接近最小值。

再來,我們不知道最佳解位置x*。我們無法求得初始誤差。我們只好將初始誤差e₀ = x* - x₀改成初始殘差r₀ = Ae₀ = b - Ax₀。整個式子乘上A,投影對象也乘上A,整個演算法過程都乘上A。

每回合將殘差投影到新建立的Aᵏb、扣除分量,使得殘差越來越小。

演算法步驟:

x₀ = random vector                 // initial position
r₀ = b - Ax₀                       // initial residual
if (‖r₀‖ < ε) return x₀            // ‖r₀‖ = sqrt(❮r₀, r₀❯)
for (k = 0 ... enough times)       // loop enough times
{
    qk = A rk                      // create Krylov basis
    αk = ❮rk, qk❯ / ❮qk, qk❯       // project to Krylov basis
    rk+1 = rk - αk qk              // update residual
    xk+1 = xk + αk rk              // update position
    if (‖rk+1‖ < ε) return xk+1    // ‖rk+1‖ = sqrt(❮rk+1, rk+1❯)
}
return xk+1

順帶一提,共軛方向法可以視作最小殘差法的改良版本。雙共軛梯度法可以取代最小殘差法。