Multivariate Function

Multivariate Function(Nonlinear Function)

「多變量函數」。輸入輸出推廣成向量。

多個函數。共用輸入。

{ f₁(x, y) = x + y
{ f₂(x, y) = 2x² + 1
{ f₃(x, y) = 1 / y

數學家重新解讀為:一個函數,輸入與輸出推廣成多個數值。微觀是多個函數,宏觀是一個函數。

F(x, y) = (x + y, 2x² + 1, 1 / y)

多個數值甚至組成向量、矩陣。多個數值視作一個元件,佯裝成普通的函數。「多個函數」佯裝成「一個多變量函數」。寫來簡便,讀來卻傷腦筋,初學者可能要花時間習慣一下。

             [  x + y  ]
F( [ x ] ) = [ 2x² + 1 ]
   [ y ]     [  1 / y  ]

   F(x)    =      y

方便起見,我們重新設定變數名稱,依序標上編號。數學家從1開始,計算學家從0開始。因為這是數學式子,所以從1開始。

              [ x₁ + x₂  ]
F( [ x₁ ] ) = [ 2x₁² + 1 ]
   [ x₂ ]     [ 1 / x₂   ]

   F(x)     =      y

Multivariate Root Finding

Multivariate Root Finding(Nonlinear Equation)

「多變量求根」。多變量函數求根。

輸出變成多個數值,衍生一個問題:如何求斜率?

多變量函數微分

函數對各個變數微分。得到各個導數。

                              ∂f               ∂f         ∂f     
f(x) = (2x₁² + x₁x₂ + x₃²)    ――― = 4x₁ + x₂   ――― = x₁   ――― = 2x₃
                              ∂x₁              ∂x₂        ∂x₃    

函數對向量微分。導數們從上到下依序併成向量。

                                  [ x₁ ]
f(x) = (2x₁² + x₁x₂ + x₃²)    x = [ x₂ ]
                                  [ x₃ ]

        ∂         [ ∂f/∂x₁ ]   [ 4x₁ + x₂ ]
f′(x) = ―― f(x) = [ ∂f/∂x₂ ] = [ x₁       ]
        ∂x        [ ∂f/∂x₃ ]   [ 2x₃      ]

多個函數,從上到下依序併成一個向量,成為多變量函數。

多變量函數對向量微分。導數們從左到右依序併成矩陣。

       [ f₁(x) ]   [ 2x₁² + x₁x₂ + x₃²  ]
       [ f₂(x) ]   [ x₁x₂ + x₂x₃ + x₃x₁ ]
F(x) = [ f₃(x) ] = [ x₁² + x₂² + x₃²    ]
       [ f₄(x) ]   [ x₁ + x₂ + x₃       ]
       [ f₅(x) ]   [ x₁ - x₂            ]

        ∂         [   |      |      |      |      |    ]
F′(x) = ―― F(x) = [ f₁′(x) f₂′(x) f₃′(x) f₄′(x) f₅′(x) ]
        ∂x        [   |      |      |      |      |    ]

                  [ 4x₁+x₂ x₂+x₃ 2x₁ 1  1 ]
                = [ x₁     x₁+x₃ 2x₂ 1 -1 ]
                  [ 2x₃    x₁+x₂ 2x₃ 1  0 ]

多變量函數微分:微分公式

微分乘法律:前微後不微,前不微後微,兩者相加。

∂                 ∂                  ∂
―― F(x)ᵀG(x) := [ ―― F(x) ] G(x) + [ ―― G(x) ] F(x)
∂x                ∂x                 ∂x

∂             ∂               ∂
―― xᵀF(x) = [ ―― x ] F(x) + [ ―― F(x) ] x = F(x) + F′(x)x
∂x            ∂x              ∂x

∂            ∂                  ∂
―― ‖F(x)‖² = ―― F(x)ᵀF(x) = 2 [ ―― F(x) ] F(x) = 2F′(x)F(x)
∂x           ∂x                 ∂x

微分連鎖律:外微,內微,從右到左連乘。

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

∂            ∂       ∂
―― ‖F(x)‖² = ―― F(x) ――――― (F(x))² = (F′(x))(2F(x)) = 2F′(x)F(x)
∂x           ∂x      ∂F(x)

多變量函數微分:一覽表

多變量函數微分的結果,等於Jacobian Matrix的轉置矩陣。

F′(x) = Jᵀ(x)

多變量函數微分、一次函數微分,兩者公式互相對應。F(x)對應Ax,Jᵀ(x)對應Aᵀ。

multivariate function     |  linear function     |  polynomial func.
——————————————————————————|——————————————————————|——————————————————
∂                         |  ∂                   |  d               
―― F(x) = Jᵀ(x)           |  ―― Ax = Aᵀ          |  ―― ax = a       
∂x                        |  ∂x                  |  dx              
                          |                      |                  
∂                         |  ∂                   |  d               
―― (F(x) + b) = Jᵀ(x)     |  ―― (Ax + b) = Aᵀ    |  ―― (ax + b) = a 
∂x                        |  ∂x                  |  dx              
                          |                      |                  
∂                         |  ∂                   |  d               
―― yᵀF(x) = Jᵀ(x)y        |  ―― yᵀAx = Aᵀy       |  ―― axy = ay     
∂x                        |  ∂x                  |  dx              
                          |                      |                  
∂                         |  ∂                   |  d               
―― xᵀF(x) = F(x) + Jᵀ(x)x |  ―― xᵀAx = (A + Aᵀ)x |  ―― ax² = 2ax    
∂x                        |  ∂x                  |  dx              
                          |                      |                  
∂                         |  ∂                   |  d               
―― ‖F(x)‖² = 2Jᵀ(x)F(x)   |  ―― ‖Ax‖² = 2AᵀAx    |  ―― (ax)² = 2a²x 
∂x                        |  ∂x                  |  dx              

多個函數的一次微分(Transpose Jacobian Matrix)

多個函數,從上到下依序併成一個向量。

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

該矩陣是「多個函數的一階導數矩陣」。通常不是方陣。

該矩陣也是Jacobian Matrix的轉置矩陣。

    [ f₁ ]
F = [ f₂ ]
    [ f₃ ]

     ∂F   [    |      |      |   ]   [ ∂f₁/∂x₁ ∂f₂/∂x₁ ∂f₃/∂x₁ ]
F′ = ―― = [ ∂f₁/∂x ∂f₂/∂x ∂f₃/∂x ] = [ ∂f₁/∂x₂ ∂f₂/∂x₂ ∂f₃/∂x₂ ]
     ∂x   [    |      |      |   ]   [ ∂f₁/∂x₃ ∂f₂/∂x₃ ∂f₃/∂x₃ ]

          [ ―――――― ∂f₁/∂x ―――――― ]ᵀ  [ ∂f₁/∂x₁ ∂f₁/∂x₂ ∂f₁/∂x₃ ]ᵀ
   = Jᵀ = [ ―――――― ∂f₂/∂x ―――――― ] = [ ∂f₂/∂x₁ ∂f₂/∂x₂ ∂f₂/∂x₃ ]
          [ ―――――― ∂f₃/∂x ―――――― ]   [ ∂f₃/∂x₁ ∂f₃/∂x₂ ∂f₃/∂x₃ ]

一個函數的二次微分(Hessian Matrix)

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

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

該矩陣是「一個函數的二階導數矩陣」。一定是對稱矩陣。

該矩陣稱作Hessian Matrix。

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

     ∂f′   [       |             |             |      ]
f″ = ――― = [ ∂(∂f/∂x₀)/∂x  ∂(∂f/∂x₁)/∂x  ∂(∂f/∂x₂)/∂x ]
     ∂x    [       |             |             |      ]

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

二次微分=原函數的Hessian Matrix=一階導數的Transpose Jacobian Matrix。由於是對稱矩陣,可以省略轉置。

f″ = Hf = Jf′ᵀ = Jf′

演算法(ℝ⇨ℝ Newton's Method)

從頭複習一遍牛頓法吧。

座標(x, f(x))的切線,與座標軸相交之處,做為下一個x。

斜率是直除以橫。得到一道等式。

高度f(xₜ),除以跨距(xₜ - xₜ₊₁),得到切線斜率f′(xₜ)。

f′(xₜ) = f(xₜ) / (xₜ - xₜ₊₁)

只有乘法,沒有除法。

f′(xₜ) (xₜ - xₜ₊₁) = f(xₜ)

移項整理,提出xₜ₊₁。得到牛頓法公式。

xₜ₊₁ = xₜ - f(xₜ) / f′(xₜ)

時間複雜度是O(T),T是遞推次數。

演算法(ℝⁿ⇨ℝ Newton's Method)

牛頓法,推廣成許多個輸入變數。

斜率推廣成梯度。切線推廣成切面。交點推廣成交線。

k[x_, y_] := Sqrt[(1.5*x*x*x)^2] + Sqrt[(2*y*y*y)^2] + 0.5*x*y - 0.05; f[x_, y_] := k[x + .4, y + .3]; x0 = .5 - .4; y0 = .65 - .3; f0 = f[x0,y0]; gx0 = N[Derivative[1,0][f][x0,y0]]; gy0 = N[Derivative[0,1][f][x0,y0]]; xyplane = Graphics3D[{Opacity[0], Polygon[{{-1,-1,0},{-1,1,0},{1,1,0},{1,-1,0}}]}, Boxed -> False (*, ViewPoint -> {0, 0, -Infinity}*) ]; func = Plot3D[f[x,y], {x, -1, 1}, {y, -1, 1}, PlotRange -> {-1, 1}, ClippingStyle -> None, Mesh -> None, PlotStyle->Opacity[0.8], RegionFunction -> Function[{x, y, z}, 0 <= z <= 1], ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]; plane = Plot3D[gx0*(x-x0)+gy0*(y-y0)+f0, {x, -1, 1}, {y, -1, 1}, Mesh -> None, PlotStyle -> Opacity[0.8], ColorFunction -> (RGBColor[0,112,192] &), RegionFunction -> Function[{x, y, z}, 0 <= z <= 1]]; root = Plot3D[0, {x, -1, 1}, {y, -1, 1}, PlotRange -> {-1, 1}, PlotStyle -> Opacity[0], Mesh -> None, RegionFunction -> Function[{x, y, z}, f[x,y] > 0]]; levelset = SliceContourPlot3D[f[x,y], z == 0, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, PlotRange -> {-1, 1}, Contours -> Table[i,{i,0,1,0.1}], ContourShading -> None, ContourStyle -> LightGray, BoundaryStyle -> LightGray]; x0bar = Graphics3D[{Black, Line[{{x0, -1, 0},{x0, y0, 0}}]}]; y0bar = Graphics3D[{Black, Line[{{-1, y0, 0},{x0, y0, 0}}]}]; f0bar = Graphics3D[{Thickness[0.01], CapForm["Butt"], Black, Line[{{x0, y0, 0},{x0, y0, f0}}]}]; f0point = Graphics3D[{Black, Sphere[{x0, y0, f0}, 0.05]}]; x1 = x0 - f0 * gx0 / ((gx0 ^ 2) + (gy0 ^ 2)); y1 = y0 - f0 * gy0 / ((gx0 ^ 2) + (gy0 ^ 2)); x1bar = Graphics3D[{Black, Line[{{x1, -1, 0},{x1, y1, 0}}]}]; y1bar = Graphics3D[{Black, Line[{{-1, y1, 0},{x1, y1, 0}}]}]; line = Graphics3D[{Black, Line[{{x0, y0, f0},{x1, y1, 0}}]}]; xs = -f0 / gx0 + x0; ys = -f0 / gy0 + y0; xsline = Graphics3D[{Black, Line[{{xs, y0, 0},{x0, y0, f0}}]}]; ysline = Graphics3D[{Black, Line[{{x0, ys, 0},{x0, y0, f0}}]}]; xspoint = Graphics3D[{Black, Sphere[{xs, y0, 0}, 0.05]}]; yspoint = Graphics3D[{Black, Sphere[{x0, ys, 0}, 0.05]}];
Show[xyplane,func,plane,f0point,f0bar,x0bar,y0bar] (* RGBColor[168,206,234] *) Show[xyplane,plane,f0point,f0bar,x0bar,y0bar,x1bar,y1bar,line] Show[xyplane,root,plane]

X方向切線與XY平面的交點、Y方向切線與XY平面的交點,兩交點連線,得到交線。

Show[xyplane,plane,f0point,f0bar,x0bar,y0bar,xsline,ysline,xspoint,yspoint] Show[xyplane,func,plane] (* ViewPoint -> {0, 0, -Infinity} *)

切面上有許多切線,可以做為前進方向;交線上有許多點,可以做為下一個x。牛頓法選擇了最短切線:垂直於交線。

Show[xyplane,plane,f0point,f0bar,x0bar,y0bar]

最短切線一定是梯度方向。等高線的垂直方向。

Show[xyplane,plane,levelset,f0point,f0bar,x0bar,y0bar,x1bar,y1bar,line] ContourPlot[f[x,y], {x, -1, 1}, {y, -1, 1}, Frame -> False, ContourShading -> None, Contours -> Table[i,{i,0,1,0.1}]]

接下來把上述概念寫成數學式子。

乘法推廣成點積。得到所有切線。

∇f(xₜ)ᵀ (xₜ - xₜ₊₁) = f(xₜ)

除法推廣成點積反運算(虛擬反矩陣)。得到最短切線。

                             ∇f(xₜ) f(xₜ)    ∇f(xₜ) f(xₜ)
(xₜ - xₜ₊₁) = ∇f(xₜ)ᵀ⁺ f(xₜ) = ―――――――――――― = ―――――――――――――
                               ‖∇f(xₜ)‖²    ∇f(xₜ)ᵀ ∇f(xₜ)

移項整理,提出xₜ₊₁。得到牛頓法公式。

                                ∇f(xₜ) f(xₜ)
xₜ₊₁ = xₜ - ∇f(xₜ)ᵀ⁺ f(xₜ) = xₜ - ―――――――――――――
                               ∇f(xₜ)ᵀ ∇f(xₜ)

時間複雜度是O(NT),N是變數數量,T是遞推次數。

演算法(ℝⁿ⇨ℝᵐ Newton's Method)

牛頓法,推廣成多變量函數。

F′(xₜ)ᵀ (xₜ - xₜ₊₁) = F(xₜ)
(xₜ - xₜ₊₁) = F′(xₜ)ᵀ⁺ F(xₜ)
xₜ₊₁ = xₜ - F′(xₜ)ᵀ⁺ F(xₜ)

函數微分改寫成Jacobian Matrix風格。

xₜ₊₁ = xₜ - J(xₜ)⁺ F(xₜ)

時間複雜度是O(N³T),N是變數數量與函數數量較大者,T是遞推次數。

演算法(ℝ⇨ℝ Secant Method)

從頭複習一遍割線法吧。

座標(x₁, f(x₁))與(x₂, f(x₂))的割線,與座標軸相交之處,做為下一個x。

牛頓法是切線斜率f′(xₜ)。割線法是割線斜率sₜ。

Newton's Method
xₜ₊₁ = xₜ - f(xₜ) / f′(xₜ)

Secant Method
{ sₜ = (f(xₜ₋₁) - f(xₜ)) / (xₜ₋₁ - xₜ)
{ xₜ₊₁ = xₜ - f(xₜ) / sₜ

只有乘法,沒有除法。

{ sₜ (xₜ₋₁ - xₜ) = f(xₜ₋₁) - f(xₜ)
{ sₜ (xₜ - xₜ₊₁) = f(xₜ)

移項整理,提出xₜ₊₁。得到割線法公式。

xₜ₊₁ = xₜ - f(xₜ) / sₜ
xₜ₊₁ = xₜ - f(xₜ) (xₜ₋₁ - xₜ) / (f(xₜ₋₁) - f(xₜ))

改寫成跨距風格。

xₜ₊₁ = xₜ - f(xₜ) Δxₜ / Δfₜ

時間複雜度是O(T),T是遞推次數。

演算法(ℝⁿ⇨ℝ Secant Method)

割線法,推廣成許多個輸入變數。

斜率推廣成梯度。割線推廣成割面。交點推廣成交線。

k[x_, y_] := Sqrt[(1.5*x*x*x)^2] + Sqrt[(2*y*y*y)^2] + 0.5*x*y - 0.05; f[x_, y_] := k[x + .4, y + .3]; x0 = .5 - .4; y0 = .65 - .3; f0 = f[x0,y0]; x1 = .425 - .4; y1 = .45 - .3; f1 = f[x1,y1]; dx = x0 - x1; dy = y0 - y1; df = f0 - f1; d = (dx ^ 2) + (dy ^ 2); gx1 = df * dx / d; gy1 = df * dy / d; g1 = (gy1 ^ 2) + (gx1 ^ 2); x2 = x1 - f1 * gx1 / g1; y2 = y1 - f1 * gy1 / g1; t[x_, y_] := gx1*(x-x1) + gy1*(y-y1) + f1; x0bar = Graphics3D[{Black, Line[{{x0, -1, 0},{x0, y0, 0}}]}]; y0bar = Graphics3D[{Black, Line[{{-1, y0, 0},{x0, y0, 0}}]}]; f0bar = Graphics3D[{Thickness[0.01], CapForm["Butt"], Black, Line[{{x0, y0, 0},{x0, y0, f0}}]}]; f0point = Graphics3D[{Black, Sphere[{x0, y0, f0}, 0.05]}]; x1bar = Graphics3D[{Black, Line[{{x1, -1, 0},{x1, y1, 0}}]}]; y1bar = Graphics3D[{Black, Line[{{-1, y1, 0},{x1, y1, 0}}]}]; f1bar = Graphics3D[{Thickness[0.01], CapForm["Butt"], Black, Line[{{x1, y1, 0},{x1, y1, f1}}]}]; f1point = Graphics3D[{Black, Sphere[{x1, y1, f1}, 0.05]}]; x2bar = Graphics3D[{Black, Line[{{x2, -1, 0},{x2, y2, 0}}]}]; y2bar = Graphics3D[{Black, Line[{{-1, y2, 0},{x2, y2, 0}}]}]; line = Graphics3D[{Black, Line[{{x0, y0, f0},{x2, y2, 0}}]}]; xyplane = Graphics3D[{Opacity[0], Polygon[{{-1,-1,0},{-1,1,0},{1,1,0},{1,-1,0}}]}, Boxed -> False (*, ViewPoint -> {0, 0, -Infinity}*) ]; func = Plot3D[f[x,y], {x, -1, 1}, {y, -1, 1}, PlotRange -> {-1, 1}, ClippingStyle -> None, Mesh -> None, PlotStyle->Opacity[0.8], RegionFunction -> Function[{x, y, z}, 0 <= z <= 1], ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]; plane = Plot3D[t[x,y], {x, -1, 1}, {y, -1, 1}, Mesh -> None, PlotStyle -> Opacity[0.8], ColorFunction -> (RGBColor[0,112,192] &), RegionFunction -> Function[{x, y, z}, 0 <= z <= 1]]; intersect = Plot3D[t[x,y]+0.02, {x, -1, 1}, {y, -1, 1}, Mesh -> None, PlotStyle -> Opacity[0.8], ColorFunction -> (RGBColor[0,112,192] &), RegionFunction -> Function[{x, y, z}, z >= f[x,y]+0.02]];
Show[xyplane,func,plane,intersect,x0bar,y0bar,f0bar,f0point,x1bar,y1bar,f1bar,f1point] (* RGBColor[168,206,234] *) Show[xyplane,plane,x0bar,y0bar,f0bar,f0point,x1bar,y1bar,f1bar,f1point,x2bar,y2bar,line]

兩點之間有許多割面。割線法選擇了最正割面:其交線垂直於兩點連線。

par = LinearSolve[{{x0, y0, 1}, {x1, y1, 1}, {x0 - 0.07, 0, 1}}, {f0, f1, 0}]; par = LinearSolve[{{x0, y0, 1}, {x1, y1, 1}, {x0 + 0.005, 0, 1}}, {f0, f1, 0}]; planebad = Plot3D[(par[[1]]*x)+(par[[2]]*y)+par[[3]], {x, -1, 1}, {y, -1, 1}, Mesh -> None, PlotStyle -> Opacity[0.8], ColorFunction -> (RGBColor[0,112,192] &), RegionFunction -> Function[{x, y, z}, 0 <= z <= 1]]; intersectbad = Plot3D[(par[[1]]*x)+(par[[2]]*y)+par[[3]]+0.02, {x, -1, 1}, {y, -1, 1}, Mesh -> None, PlotStyle -> Opacity[0.8], ColorFunction -> (RGBColor[0,112,192] &), RegionFunction -> Function[{x, y, z}, 1 >= z >= f[x,y]+0.02]]; Show[xyplane,func,planebad,intersectbad,x0bar,y0bar,f0bar,f0point,x1bar,y1bar,f1bar,f1point]

割面上有許多割線,可以做為前進方向;交線上有許多點,可以做為下一個x。割線法選擇了最短割線:垂直於交線。

Show[xyplane,plane,x0bar,y0bar,f0bar,f0point,x1bar,y1bar,f1bar,f1point]

接下來把上述概念寫成數學式子。

點積。得到所有割面。

sₜᵀ (xₜ₋₁ - xₜ) = f(xₜ₋₁) - f(xₜ)

點積反運算(虛擬反矩陣)。得到最正割面。

sₜᵀ = (f(xₜ₋₁) - f(xₜ)) (xₜ₋₁ - xₜ)⁺ = Δfₜ Δxₜ⁺

點積。得到所有割線。

sₜᵀ (xₜ - xₜ₊₁) = f(xₜ)

點積反運算(虛擬反矩陣),得到最短割線。

(xₜ - xₜ₊₁) = sₜᵀ⁺ f(xₜ) = (Δfₜ Δxₜ⁺)⁺ f(xₜ) = Δxₜ Δfₜ⁺ f(xₜ)

移項整理,提出xₜ₊₁。得到割線法公式。

xₜ₊₁ = xₜ - Δxₜ Δfₜ⁺ f(xₜ)

然而割線法只會走一直線,往往走不到根。沒人使用割線法。

演算法(ℝⁿ⇨ℝᵐ Secant Method)

割線法,推廣成多變量函數。

xₜ₊₁ = xₜ - Δxₜ ΔFₜ⁺ F(xₜ)

一樣沒人用。

演算法(Broyden's Method)

割線法改良版。不選擇最正割面,避免走一直線。

割線法寫成兩道式子,先改良第一式,再改良第二式。

{ Sₜᵀ = (F(xₜ₋₁) - F(xₜ)) (xₜ₋₁ - xₜ)⁺ = ΔFₜ Δxₜ⁺
{ xₜ₊₁ = xₜ - Sₜᵀ⁺ F(xₜ)

一、改成遞推更新割面。用前一個割面,求後一個割面。

故意製造Sₜ₋₁。移項整理,提出Sₜ。得到割面遞迴公式。

最初的割面,設定成隨機梯度(隨便亂歪),或者第一點的切面的梯度(不會太歪)。

最初的割面,故意設定成歪的;之後的割面,故意保持是歪的。得以螺旋行進至根。

Sₜᵀ Δxₜ = ΔFₜ
(Sₜ - Sₜ₋₁ + Sₜ₋₁)ᵀ Δxₜ = ΔFₜ
(Sₜ - Sₜ₋₁)ᵀ Δxₜ + Sₜ₋₁ᵀ Δxₜ = ΔFₜ
(Sₜ - Sₜ₋₁)ᵀ Δxₜ = ΔFₜ - Sₜ₋₁ᵀ Δxₜ
(Sₜ - Sₜ₋₁)ᵀ = (ΔFₜ - Sₜ₋₁ᵀ Δxₜ) Δxₜ⁺
Sₜᵀ - Sₜ₋₁ᵀ = (ΔFₜ - Sₜ₋₁ᵀ Δxₜ) Δxₜ⁺     Δxₜ⁺ = Δxₜᵀ / ‖Δxₜ‖² 
Sₜᵀ = Sₜ₋₁ᵀ + (ΔFₜ - Sₜ₋₁ᵀ Δxₜ) Δxₜ⁺         = Δxₜᵀ / (Δxₜᵀ Δxₜ)

第一式改良完畢,得到計算公式。

{ Sₜᵀ = Sₜ₋₁ᵀ + (ΔFₜ - Sₜ₋₁ᵀ Δxₜ) Δxₜ⁺
{ xₜ₊₁ = xₜ - Sₜᵀ⁺ F(xₜ)

二、改成直接紀錄虛擬反矩陣。避免不斷計算虛擬反矩陣。

令Cₜᵀ = Sₜᵀ⁺。

Sₜᵀ Δxₜ = ΔFₜ
Cₜᵀ ΔFₜ = Δxₜ
(Cₜ - Cₜ₋₁ + Cₜ₋₁)ᵀ ΔFₜ = Δxₜ
(Cₜ - Cₜ₋₁)ᵀ ΔFₜ + Cₜ₋₁ᵀ ΔFₜ = Δxₜ
(Cₜ - Cₜ₋₁)ᵀ ΔFₜ = Δxₜ - Cₜ₋₁ᵀ ΔFₜ
(Cₜ - Cₜ₋₁)ᵀ = (Δxₜ - Cₜ₋₁ᵀ ΔFₜ) ΔFₜ⁺
Cₜᵀ - Cₜ₋₁ᵀ = (Δxₜ - Cₜ₋₁ᵀ ΔFₜ) ΔFₜ⁺     ΔFₜ⁺ = ΔFₜᵀ / ‖ΔFₜ‖²
Cₜᵀ = Cₜ₋₁ᵀ + (Δxₜ - Cₜ₋₁ᵀ ΔFₜ) ΔFₜ⁺         = ΔFₜᵀ / (ΔFₜᵀ ΔFₜ)

第二式改良完畢,得到計算公式。

{ Cₜᵀ = Cₜ₋₁ᵀ + (Δxₜ - Cₜ₋₁ᵀ ΔFₜ) ΔFₜ⁺
{ xₜ₊₁ = xₜ - Cₜᵀ F(xₜ)

時間複雜度是O(NMT),N是變數數量,M是函數數量,T是遞推次數。

演算法(Fixed Point Iteration)

不動點遞推法,推廣成多變量函數。

函數和變數必須一樣多個。

xₜ₊₁ = xₜ + λ F(xₜ)

Multivariate Optimization

Multivariate Optimization(Nonlinear Least Squares)

「多變量最佳化」。多變量函數最佳化。

輸出變成多個數值,衍生兩個問題:如何比大小?如何求斜率?

多變量函數長度

比大小有兩種方式:

Multiobjective Optimization 輸出數值們的加權平均值(非長度函數)
Multivariate Optimization   輸出數值們的平方和(長度平方函數)

本文介紹平方和。目標函數是輸出數值們的平方和‖F(x)‖²。

       [ f₁(x) ]
F(x) = [ f₂(x) ]
       [ f₃(x) ]

‖F(x)‖² = F(x)ᵀF(x) = f₁(x)² + f₂(x)² + f₃(x)²

演算法(Gauss-Newton Algorithm)

牛頓法求極值。

一次微分等於零的地方是極值、鞍點。事先微分,然後求根。

xₜ₊₁ = xₜ - f″(xₜ)ᵀ⁺ f′(xₜ)

目標函數是‖F(x)‖²。省略二階導數,簡化計算過程。

xₜ₊₁ = xₜ - [ ∂²/∂x² ‖F(x)‖² ]ᵀ⁺ [ ∂/∂x ‖F(x)‖² ]     求值就不寫了
xₜ₊₁ = xₜ - [ 2 F″(xₜ)ᵀ F(xₜ) + 2 F′(xₜ) F′(xₜ)ᵀ]ᵀ⁺ [ 2 F′(xₜ) F(xₜ) ]
xₜ₊₁ ≈ xₜ - [ 2 F′(xₜ) F′(xₜ)ᵀ ]ᵀ⁺ [ 2 F′(xₜ) F(xₜ) ]   省略二階導數
xₜ₊₁ ≈ xₜ - [ F′(xₜ) F′(xₜ)ᵀ ]ᵀ⁺ F′(xₜ) F(xₜ)           分子分母同除2
xₜ₊₁ ≈ xₜ - [ F′(xₜ) F′(xₜ)ᵀ ]⁺ F′(xₜ) F(xₜ)            對稱矩陣,省略轉置

函數微分改寫成Jacobian Matrix風格。

xₜ₊₁ ≈ xₜ - [ Jᵀ(xₜ) Jᵀ(xₜ)ᵀ ]⁺ Jᵀ(xₜ) F(xₜ)
xₜ₊₁ ≈ xₜ - [ Jᵀ(xₜ) J(xₜ) ]⁺ Jᵀ(xₜ) F(xₜ)     ---> xₜ₊₁ = xₜ - (JᵀJ)⁺JᵀF

碰巧是Normal Equation,精簡成虛擬反矩陣。

xₜ₊₁ ≈ xₜ - J(xₜ)⁺ F(xₜ)     ---> xₜ₊₁ = xₜ - J⁺F

多變量函數,牛頓法求根、牛頓法求極值,公式一模一樣。

演算法(Gradient Descent)

不動點遞推法求極值。

xₜ₊₁ = xₜ - λ f′(xₜ)

目標函數是‖F(x)‖²。

xₜ₊₁ = xₜ - λ [ ∂/∂x ‖F(x)‖² ]
xₜ₊₁ = xₜ - λ [ 2 F′(xₜ) F(xₜ) ]
xₜ₊₁ = xₜ - λ F′(xₜ) F(xₜ)           倍率2併入步伐大小λ

函數微分改寫成Jacobian Matrix風格。

xₜ₊₁ = xₜ - λ Jᵀ(xₜ) F(xₜ)      ---> xₜ₊₁ = xₜ - λJᵀF

演算法(Levenberg-Marquardt Algorithm)

牛頓法求極值:前進方向因地制宜,攻頂路線短。步伐大小不可調整(二階導數的倒數,曲率半徑越大、路線越直、步伐越大),最初只能亦步亦趨。

不動點遞推法求極值:前進方向不可調整(梯度方向),攻頂路線長。步伐大小可以調整,最初就能一步登天。

Gauss-Newton Algorithm與Gradient Descent的前進方向,取加權平均。隨時視情況調整權重,兩害相權取其輕。

xₜ₊₁ = xₜ - (J⁺ + λJᵀ)F
            ^^   ^^^
        Newton   fixpoint

有人將Gradient Descent做梯度正規化。

JᵀJ的對角線元素,是梯度的每個元素的平方。diag(JᵀJ)將之併成一個向量。diag(JᵀJ)⁺形成倒數,即是梯度正規化──然而他忘了開根號。大家將錯就錯。

xₜ₊₁ = xₜ - (J⁺ + diag(JᵀJ)⁺ λJᵀ)F
            ^^   ^^^^^^^^^^ ^^^
        Newton   normalize  fixpoint

有人從Normal Equation那道式子下手,寫成另外一種風格。

xₜ₊₁ = xₜ - [JᵀJ + 1/λ diag(JᵀJ)]⁺ JᵀF
            ^^^   ^^^^^^^^^^^^^
         Newton   normalized fixpoint

牛頓法公式一覽表

虛擬反矩陣改寫成矩陣除法風格。

ℝ⇨ℝ   求根  xₜ₊₁ = xₜ -  f(xₜ) / f′(xₜ)
ℝ⇨ℝ   求極值 xₜ₊₁ = xₜ - f′(xₜ) / f″(xₜ)
ℝⁿ⇨ℝ  求根  xₜ₊₁ = xₜ -  f(xₜ) / f′(xₜ)ᵀ = xₜ - J⁺f
ℝⁿ⇨ℝ  求極值 xₜ₊₁ = xₜ - f′(xₜ) / f″(xₜ)ᵀ = xₜ - H⁺f′
ℝⁿ⇨ℝᵐ 求根  xₜ₊₁ = xₜ -  F(xₜ) / F′(xₜ)ᵀ = xₜ - J⁺F
ℝⁿ⇨ℝᵐ 求極值 xₜ₊₁ = xₜ - E′(xₜ) / E″(xₜ)ᵀ ≈ xₜ - J⁺F , E(x) = ‖F(x)‖²