Visual Optimization in two dimensions

Written by Adam B. Levy, April 2009

This notebook gives a template for visually observing the progress of different numerical methods for optimizing a function f(x,y) of two-variables, possibly subject to a constraint g(x,y)≥0.

Intializing the Guess and Objective

Off[General::spell]
Off[General::spell1]
Clear[f]
nb=SelectedNotebook[ ];
inx=.4;
iny=.7;
fi=x^3-x y^2+y^3-y+1;
gi=0;
(*The next commands build partial derivatives that can be used in the iteration function (and an exact line search) while accepting current guesses in the argument later*)
fy[x_,y_]=D[f[x,y],y];
fx[x_,y_]=D[f[x,y],x];
fxx[x_,y_]=D[fx[x,y],x];
fxy[x_,y_]=D[fx[x,y],y];
fyx[x_,y_]=D[fy[x,y],x];
fyy[x_,y_]=D[fy[x,y],y];
PasteButton[Defer[fx[x,y]]];
PasteButton[Defer[fy[x,y]]];
PasteButton[Defer[fxx[x,y]]];
PasteButton[Defer[fxy[x,y]]];
PasteButton[Defer[fyx[x,y]]];
PasteButton[Defer[fyy[x,y]]];
(*The next command finds the exact line search step of f from {xo,yo} in the direction {d1,d2}*)
ExactStep[xo_,yo_,{d1_,d2_}]:=Module[{a,x,y,min},min=Minimize[f[xo+a d1,yo+a d2],a];a/.min[[2]]]
(*The next command finds the solution to the SQP quadratic program with multiplier lam at (xo,yo)*)SQP[xo_,yo_,lam_]:=Module[{L,x,y,min,xn,yn},L[x_,y_]=f[x,y]-lam g[x,y];min=Minimize[{L[xo,yo]+(D[L[x,y],x]/.{x,y}->{xo,yo})(x-xo)+(D[L[x,y],y]/.{x,y}->{xo,yo})(y-yo)+1/2{x-xo,y-yo}.{{(D[L[x,y],x,x]/.{x,y}->{xo,yo}),(D[L[x,y],x,y]/.{x,y}->{xo,yo})},{(D[L[x,y],y,x]/.{x,y}->{xo,yo}),(D[L[x,y],y,y]/.{x,y}->{xo,yo})}}.{x-xo,y-yo},g[xo,yo]+(D[g[x,y],x]/.{x,y}->{xo,yo})(x-xo)+(D[g[x,y],x]/.{x,y}->{xo,yo})(y-yo)≥0},{x,y}];{Replace[x,min[[2]]],Replace[y,min[[2]]]}]
(*The next command iterates the iteration function once and updates the list of guesses*)
IterateOpt[{xo_,yo_}]:=Module[{xn,yn},{xn,yn}=Evaluate[G[xo,yo]];guess=AppendTo[guess,N[{xn,yn}]];]
(*The next command draws the contour plot of the objective with the guesslist and constraint curve--if there is one--overlayed*)
VisualOpt[list_,f_,g_]:=Module[{xmax,xmin,ymax,ymin},xmax=Max[list[[All,1]]];xmin=Min[list[[All,1]]];ymax=Max[list[[All,2]]];ymin=Min[list[[All,2]]];
Grid[{{Grid[{{"newest 20 guesses"},{If[Length[list]≤19,TableForm[list],TableForm[Take[list,-20]]]}},Alignment->{Center,Center}],Show[ContourPlot[f,{x,xmin-If[xmax==xmin,1,(xmax-xmin)]/10,xmax+If[xmax==xmin,1,(xmax-xmin)]/10},{y,ymin-If[ymax==ymin,1,(ymax-ymin)]/10,ymax+If[ymax==ymin,1,(ymax-ymin)]/10},Contours->20,PlotPoints->25],Graphics[{Blue,Line[list]}],Graphics[{Blue,PointSize[Medium],Point[list]}],Graphics[{Blue,PointSize[Large],Point[First[list]]}],Graphics[{Red,PointSize[Large],Point[Last[list]]}],ContourPlot[g,{x,xmin-If[xmax==xmin,1,(xmax-xmin)]/10,xmax+If[xmax==xmin,1,(xmax-xmin)]/10},{y,ymin-If[ymax==ymin,1,(ymax-ymin)]/10,ymax+If[ymax==ymin,1,(ymax-ymin)]/10},Contours->{0},PlotPoints->25,ContourShading->False,ContourStyle->{Thick,Blue}],Frame->True,FrameLabel->{Style["Newest guess is "Style[Last[list],Red],Medium],,Style["First guess is "Style[First[list],Blue],Medium]},ImageSize->Medium]}},ItemSize->{{15,30}},Alignment->{Center,Top},Frame->All]]
Button["Step",IterateOpt[Last[guess]];Print[VisualOpt[guess,fi,gi]]];
(*The next command makes a pop-up box to intiialize the method*)
Initialize:=Module[{},DialogInput[{Grid[{{TextCell["First Guess"]},{Grid[{{TextCell["x ="],InputField[Dynamic[inx],Number,ImageSize->30],TextCell["y ="],InputField[Dynamic[iny],Number,ImageSize->30]}}]}}],Grid[{{TextCell["Objective function"]},{TextCell["f[x,y]"],InputField[Dynamic[fi]]}}],Grid[{{TextCell["Constraint function"]},{TextCell["g[x,y]"],InputField[Dynamic[gi]],TextCell["≥ 0"]}}],
DefaultButton["Done",DialogReturn[guess={{inx,iny}};f[x_,y_]=fi;g[x_,y_]=gi;]]}]]

To start the numerical method at a particular guess with a particular objective function, click your cursor after the following "Initialize" command, hit the ENTER button on your keyboard (not the carriage return button), click "Yes" to the initialization question, fill in the appropriate boxes that pop up, and click "Done".

Initialize

Building the Method

Iteration Function

The numerical method is determined by an iteration function G(x,y) giving the right side of the update formula for new guesses in terms of old:

[({{(Subscript[x, new])}, {(Subscript[y, new])}})] = G((Subscript[x, old]),(Subscript[y, old]))

For example, simply assigning the old variables as new would involve the iteration function

G(x,y)={x,y}

since Mathematica represents column vectors [({{x}, {y}})] as lists {x,y}.

ENTER your iteration function after filling in the □ below using the following pasting buttons

({{fx[x,y], fy[x,y], fxx[x,y], fxy[x,y], fyx[x,y], fyy[x,y]}})

G[x_,y_]:=\[Placeholder]

Open this cell for Mathematica hints on vector/matrix multiplication

A vector (({{x}, {y}})) can be represented in Mathematica as: {x,y}
A matrix (({{a, b}, {c, d}})) can be represented in Mathematica as: {{a,b},{c,d}}
An inverse matrix (((({{a, b}, {c, d}})))^(-1)) can be represented in Mathematica as:  Inverse[{{a,b},{c,d}}]
A matrix-vector multiplication (({{a, b}, {c, d}})) (({{x}, {y}})) can be represented in Mathematica as:  {{a,b},{c,d}}.{x,y}

Open this cell for a command to use for an exact line search

ExactStep[x,y,{d1,d2}] returns the step-factor α from the exact line search minimizing f(x,y) in the step-direction [({{d1}, {d2}})] from the base point (x,y)

Open this cell for a command to solve the QP in SQP

SQP[xold,yold,λold] returns the solution (x,y) to the quadratic program from the SQP method at the point (xold,yold) with multiplier λold.

Computing Guesses

Press the step button below to update the new guess from the old.  Repeat as often as you like.

"Step"

Open this cell for details on how to extract old guesses from the full list

If at any time, you need to compute something else involving the guesses, you can extract them from the list guess in which they are stored.  For example, the x and y coordinates respectively of the newest guess are obtained by entering

Last[guess][[1]]
Last[guess][[2]]

The function Last retrieves the last element of the list (of pairs), and the [[1]] or [[2]] identifies the first or second component of the pair.