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ₜ⁺
Sₜᵀ = Sₜ₋₁ᵀ + (ΔFₜ - Sₜ₋₁ᵀ Δxₜ) Δxₜ⁺     Δxₜ⁺ = Δ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ₜ⁺
Cₜᵀ = Cₜ₋₁ᵀ + (Δxₜ - Cₜ₋₁ᵀ ΔFₜ) ΔFₜ⁺     ΔFₜ⁺ = Δ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)‖²

Level Set(Under Construction!)

Level Set

多變量函數是箭矢圖。據我所知似乎沒有等高線的概念。

容易計算的多變量函數類型

遞增函數、凸函數、平緩函數都必須重新定義。我找不到相關資料,事情相當棘手。

《Learning Submodular Functions》
Submodularity is in many ways similar to concavity of functions
defined on Rn. For example, concavity of differentiable functions
is equivalent to gradient monotonicity, and submodularity is equivalent
to monotonicity of the marginal values:

• Non-negative if f(S) ≥ 0 for all S.
• Monotone (or non-decreasing) if f(S) ≤ f(T) for all S ⊆ T.
• Submodular if f(A) + f(B) ≥ f(A ∪ B) + f(A ∩ B) ∀A,B ⊆ [n].
  or equivalently f(A ∪ {i}) − f(A) ≥ f(B ∪ {i}) − f(B)
  for all A ⊆ B ⊆ [n] and i ∉ B.
• Subadditive if f(S) + f(T) ≥ f(S ∪ T) for all S,T ⊆ [n].
• L-Lipschitz if |f(S + x) − f(S)| ≤ L for all S ⊆ [n] and x ∈ [n].

Increasing Function ⮕ Nonnegative Divergence
【尚待確認】

                           f(b) - f(a)
function f is monotone iff ――――――――――― >= 0   斜率大於零,分子分母同號
                              b - a
                                          ∂   ∂
multivariate function f is monotone iff ( ――, ―― ) dot (f(b), f(a)) >= 0
                                          ∂b  ∂a
點積大於零,分子分母同號(散度大於零,被定義成向量變大。)

一維箭矢圖有個特性:遞增函數(處處斜率大於零),箭矢自零散開。遞減函數(處處斜率小於零),箭矢聚集至零。

多維箭矢圖則是:處處散度大於零,箭矢自零散開。處處散度小於零,箭矢聚集至零。

散開的地方稱作源,聚集的地方稱作匯。當函數起起伏伏、而且很多地方是零,就有許多源匯。正負無限大的地方也是源匯。

儘管源匯很吸睛,然並卵。

Increasing Function ⮕ Positive Semidefinite Matrix

                           f(b) - f(a)
function f is monotone iff ――――――――――― >= 0   斜率大於零,分子分母同號
                              b - a

multivariate function f is monotone iff (b-a) dot (f(b)-f(a)) >= 0
                                                  ^^^^^^^^^^^
                                          f(b-a) when f is linear
點積大於零,分子分母同號(轉彎幅度少於±90˚,被定義成向量變大。)

遞增數列,就是數字越來越大。遞增陣列,就是數字往右往下越大、往左往上越小。

遞增函數,就是數字越來越大。遞增二元函數,就是數字往右上的各種方向都越來越大。

輸入輸出推廣成向量,則用向量長度比大小。然而向量長度沒有負數,沒辦法越來越小。因此有人改用向量內積比大小:正數(向量轉彎幅度少於±90˚)代表更大,負數代表更小。意義扭曲了。

【輸入向量和輸出向量不一樣長,沒辦法點積,該怎麼辦?】

多變量函數是遞增函數:處處梯度(一次微分、雅可比矩陣的轉置)是半正定矩陣。處處斜率是正數或零。

線性函數Ax是遞增函數:處處梯度都是Aᵀ,Aᵀ是半正定矩陣。

嚴格遞增則是正定矩陣。特徵值都是正數,沒有零。

註:「遞增矩陣」與「單調矩陣」這些詞彙已經被古人賦予其他意義,因此大家直接稱呼為正半定矩陣。

vector field is monotone iff jacobian matrix is positive semidefinite everywhere.
linear multivariate function Ax is monotone iff A is positive semidefinite.

Lipschitz Function ⮕ Condition Number
【尚待確認】

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

Lipschitz Function ⮕ Diagonally Dominant Matrix
【尚待確認】

多變量函數的遞增函數定義成一次微分是半正定矩陣。

多變量函數的平緩函數定義成對角優勢矩陣。

monotone: (f(b) - f(a)) / (b - a) >= 0
Lipschitz: -1 <= (f(b) - f(a)) / (b - a) <= +1
           |f(a) + f(b)| <= |a + b|  其中一個變數變號
positive Lipschitz => monotone
diagonally dominant matrix = 每個橫條都是Lipschitz function
(對於橫條的非主對角線元素來說。)
(例如Gauss–Seidel,更新x,方程式5x + y = 5是陡峭的,很近就撞到,所以收斂。
(更新式x = 1 - 0.2y,對y而言,斜率-1 ≤ -0.2 ≤ 1是平緩的。)
positive diagonally dominant => positive definite => monotone
https://books.google.com.tw/books?id=KVfSBwAAQBAJ&pg=PA105