Optimization Methods

Steepest Descent

Function of One Variable

Search for the minimum of function f  by doing steps proportional to (∂_x) f in the direction, opposite to (∂_x) f

\(SteepestDescent[f_, 
        x0_, \ η_] := 
Module[
{x1 = 
            x0, n = 0, δf, 
          ptList = {x0}}, 
δf[
            x_] = ∂\_x f[
              x]; 

While[
          Abs[δf[x1]] > 0.000001, 
x1 = 
            x1 - η\ δf[
                  x1]; 
\(n++\); 
ptList = Append[ptList, x1]; 
If[
            n == 
              1000, 
Print["Didn't converge nicely"]; 
Break[]
];
]; 
Print["Converged to ", x1, " in ", 
          n, " iterations."]; 
{x1, 
          ptList}
]\)

f[x_] := (((x - 3))^2) - 4

Tmp = SteepestDescent[f, 10, 0.9]

PlotConvergence[f, Tmp[[2]]]

Tmp = SteepestDescent[f, 1, 0.3]

PlotConvergence[f, Tmp[[2]]]

SteepestDescent[f, 1, 0.5]

SteepestDescent[f, 1, 1]

SteepestDescent[f, 1, 1.1]

Function of Several Variables

Here we just use the gradient ∇finstead of the derivative

<< Calculus`VectorAnalysis`

\(SteepestDescent3D[f_, 
        x0_, \ η_] := 
Module[
{x1 = 
            x0, n = 0, δf}, 
δf[{x_, y_, 
              z_}] = Grad[f[x, y, z], 
            Cartesian[x, y, z]]; 

While[
          Max[Abs[δf[x1]]] > 0.0001, 
x1 = 
            x1 - η\ δf[
                  x1]; 
\(n++\); 
If[
            n == 
              10000, 
Print["Didn't converge nicely"]; 
Break[]
];
]; 
Print["Converged to ", x1, " in ", 
          n, " iterations."]; 
x1
]\)

f3[x_, y_, 
        z_] := (((x - 2))^2) + (((y - 3))^2) + (((z + 8))^4)

SteepestDescent3D[f3, {0, 0, 0}, 0.01]

SteepestDescent3D[f3, {0, 0, 0}, 0.001]

SteepestDescent3D[f3, {1, 1, -7.9}, 0.05]

Newton's Method

Function of One Variable

In order to find minimum of f, solve the equation ∂f = 0 using the approximation:
 f(x) ≈ 
              f((Subscript[x, 0])) + ∂f((Subscript[x, 0])) (x - (Subscript[x, 0]))+(1/2) (∂^2) f((Subscript[x, 0])) (((x - (Subscript[x, 0])))^2)  
 							⇒  ∂f(
            x) ≈ ∂f(
              (Subscript[x, 0])) + (∂^2) f((Subscript[x, 0])) (x - (Subscript[x, 0]))
Therefore.
∂f(x) = 0  ⇒  ∂f(
              (Subscript[x, 0])) + (∂^2) f((Subscript[x, 0])) (x - 
                (Subscript[x, 0]))  ≈  0  
							⇒  \(TraditionalForm\`Δ\ x\  := \ 
        x - x\_0\  ≈ \(-∂\(f(x\_0)\)\)\/∂\^2\( f(x\_0)\)\)

\(NewtonsMethod[f_, 
        x0_] := 
Module[
{x1 = x0, 
          n = 0, δf, δδf, 
          ptList = {x0}}, 
δf[
            x_] = ∂\_x f[x]; 
δδf[
            x_] = ∂\_x δf[
              x]; 

While[
          Abs[δf[x1]] > 0.000001, 
x1 = 
            x1 - N[δf[x1]\/δδf[x1]]; 
\(n++\); 
ptList = 
            Append[ptList, x1]; 
If[
            n == 
              1000, 
Print["Didn't converge nicely"]; 
Break[]
];
]; 
Print["Converged to ", x1, " in ", 
          n, " iterations."]; 
{x1, 
          ptList}
]\)

f[x_] := (((x - 100))^4)

Tmp = NewtonsMethod[f, 30]

PlotConvergence[f, Tmp[[2]]]

Function of Several Variables

The case is similar to the above, only we use the gradient ∇ finstead of the derivative, and matrix of second derivatives (kinda "gradient of a gradient") instead of the second derivative.

\(NewtonsMethod3D[f_, 
        x0_] := 
Module[
{x1 = x0, 
          n = 0, δf, δδf}, 
δf[{x_, y_, z_}] = 
          Grad[f[x, y, z], 
            Cartesian[x, y, z]]; 
δδf[{x_, 
              y_, z_}] = 
          Grad[δf[{x, y, z}], 
            Cartesian[x, y, z]]; 

While[
          Max[Abs[δf[x1]]] > 0.000001, 
x1 = 
            x1 - N[Inverse[δδf[x1]] . δf[
                    x1]]; 
\(n++\); 
If[
            n == 
              1000, 
Print["Didn't converge nicely"]; 
Break[]
];
]; 
Print["Converged to ", x1, " in ", 
          n, " iterations."]; 
x1
]\)

f3[x_, y_, 
        z_] := (((x + 6))^2) + (((y + z - 3))^2) + (((y + 9))^3)

NewtonsMethod3D[f3, {1, 2, 3}]

NewtonsMethod3D[f3, {1000, 1000, 1000}]

Gauss-Newton's Method

Suppose there is a vector-valued function e(x, y, z) ⟶ ((Subscript[e, 1]), (Subscript[e, 2]), 
          (Subscript[e, 3])) and we need to minimize the value of E = (1/2) (((Subscript[e, 1])^2) + ((Subscript[e, 2])^2) + ((Subscript[e, 3])^2)).
Linearizing e we get that for each i = 1, 2, 3
 (Subscript[e, i])(x, y, z)≈(Subscript[e, i])((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])) +  ∇(Subscript[e, i])((Subscript[x, 0]), (Subscript[y, 0]), 
          (Subscript[z, 0]))· (((x - (Subscript[x, 0]), y - (Subscript[y, 0]), z - (Subscript[z, 0])))^T)
which in vector form is:
e(x, y, z) ≈ 
        e((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])) + J·Δ(x, y, z)    where J = ∇e((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])) is the Jacobian matrix of e and Δ(x, y, z) = (x - (Subscript[x, 0]), y - (Subscript[y, 0]), 
          z - (Subscript[z, 0]))

To minimize E we solve ∇E  =  0 using the obtained approximation:
E=(1/2)[[
                  e(x, y, z)]]^2≈(1/2)[[
                    e((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0]))]]^2+(e^T)((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0]))·
                J·Δ(x, y, z)+(1/2)((Δ(x, y, z))^T)·
                  (J^T)·J·Δ(x, y, z)
∇E = 
        (J^T)·e((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])) + 
          (J^T)·J·Δ(x, y, z)

Expressing (x, y, z) we get the step of the algorithm:
\(TraditionalForm\`\((x, y, 
          z)\) = \(\((x\_0, y\_0, 
              z\_0)\) - \((J\^T·J)\)\^\(-1\)·J\^T·\(e(x\_0, y\_0, z\_0)\) = \((x\_0, y\_0, z\_0)\) - 
            J\^\(-1\)·\(e(x\_0, y\_0, z\_0)\)\)\)

\(GaussNewtonsMethod[f_, 
        x0_] := 
Module[
{x1 = x0, 
          n = 0, δf}, 
δf[{x_, y_, z_}] = 
          Grad[f[{x, y, z}], 
            Cartesian[x, y, 
              z]]; 

While[\((f[x1] . 
                f[x1])\) > 0.000001, 
x1 = 
            x1 - N[Inverse[δf[x1]] . 
                  f[x1]]; 
\(n++\); 
If[
            n == 
              1000, 
Print["Didn't converge nicely"]; 
Break[]
];
]; 
Print["Converged to ", x1, " in ", 
          n, " iterations."]; 
x1
]\)

g[{x_, y_, z_}] := {(x^2) - (y^2), (((x + z))^2) - 4 + y, 
        (z^4)}

GaussNewtonsMethod[g, {0, 0, 0}]

Tmp = GaussNewtonsMethod[g, {1, 1, 1}]

g[Tmp] . g[Tmp]

Tmp = GaussNewtonsMethod[g, {1000, -1000, 1000}]

g[Tmp] . g[Tmp]

PlotConvergence[f_, 
        list_] := 
Show[
Plot[
          f[x], {x, Min[list] - 0.2, Max[list] + 0.2}, 
          DisplayFunction -> Identity], 
ListPlot[
          Map[({#, f[#]}) &, list], PlotJoined -> True,  
          PlotStyle -> {Thickness[0.01]}, 
          DisplayFunction -> 
            Identity], 
PlotRange -> All, 
        DisplayFunction -> $DisplayFunction]

© Konstantin Tretyakov