Direction Methods.

This notebook has several demonstrations of unconstrained optimization methods which work along the following basic lines:

Given a function f:(ℝ^n)-> ℝ and a starting point (Subscript[x, 0]):

	1. Choose a direction (Subscript[v, i]) (depending on (Subscript[x, i]) and f).
	2. Find a (local) minimum of the function (Subscript[g, i])(t) = f((Subscript[x, i])+ t (Subscript[v, i])).
	3. Set (Subscript[x, i+1]) to be the location of this local min.
	
and iterate.

To make this work, we need a line search method; we’ll just use Mathematica’s FindMinimum. Note that we need to specify that we want to find a minimum in the positive direction (Subscript[v, i]), so we need t > 0.

FindMinimum[{Cos[t],t>0},{t}]

Now we can define a general direction search method:

DirectionStep[{f_,Dir_,x_}] :=Module[{v,t},
v = Dir[x];
t = FindMinimum[{f[x + t v],t > 0},{t}][[2,1,2]];
{f,Dir,x+t v}
];

DirectionMethod[f_,Dir_,x0_,n_] := NestList[DirectionStep,{f,Dir,x0},n];

Test Functions.

We’ll need some test functions to try out. First, we give a quadratic function in the special form where we can see the location of the minimum:

quadF[{x_,y_}] := ({x,y} - {2,3}).((({{12, 6}, {6, 5}})).({x,y} - {2,3}));

quadPlot = ContourPlot[Log[quadF[{x,y}]+40],{x,-20,20},{y,-20,20},Contours->20]

We also define a more complicated function with lots of peaks:

PeakLocs = {{1.26,2.22},{3.27,-3.91},{-1.54,1.67},{1.28,0.80},{-1.58,-3.72},{-1.45,-2.88},{0.95,-0.716},{1.31,0.689},{4.50,3.39},{-2.36,4.56},{1.68,1.63},{-3.44,-0.80},{1.22,-0.49},{-3.97,-2.99},{-0.42,-3.00},{-1.38,-4.91},{-3.96,1.54},{-4.38,-0.04},{-3.13,3.85},{1.91,1.95},{0.35,4.32},{-2.32,1.69},{-3.83,0.83},{-4.67,-1.19},{1.87,-1.85},{4.52,-3.29},{-1.34,-3.742},{-4.98,-3.20},{2.17,-1.83},{0.79,4.76}};

PeakHeights ={-0.257,0.12,-0.34,1.44,1.57,1.437,-0.146,0.475,0.609,1.58,-0.333,0.286,0.934,-0.0627,1.799,1.68,1.21,1.47,0.142,1.600,1.824,0.210,0.027,1.99,-0.379,0.087,0.825,1.81,1.07,0.916};

Peak[height_,{a_,b_},{x_,y_}] := N[height*Exp[-(x -a)^2 - (y -b)^2]];
peaks[{x_,y_}] := Sum[Peak[PeakHeights[[i]],PeakLocs[[i]],{x,y}],{i,1,30}] + 0.05*x^2 + 0.05*y^2;

peakPlot = ContourPlot[peaks[{x,y}],{x,-8,8},{y,-8,8},Contours->25]

Coordinate Search.

The easiest method of choosing directions is just to alternate between the directions {1,0}, {0,1}, {-1,0}, and {0,-1}.

cdcalls = 0;
CoordinateDir[x_] :=First[RotateLeft[{{1,0},{0,1},{-1,0},{0,-1}},cdcalls++]];

We can take a look at how this does over the course of 100 steps:

CoordinateSearchQuad = Last /@ DirectionMethod[quadF,CoordinateDir,{10,-17},100];

Show[quadPlot,Graphics[{Red,Line[CoordinateSearchQuad]}]]

In table form, we can see the relatively unimpressive performance of this method even more clearly:

CoordinateSearchQuad // TableForm

The problem, of course, is how we pick the directions.

CoordinateSearchPeaks = Last /@ DirectionMethod[peaks,CoordinateDir,{-4.2,-5.3},15];

Show[peakPlot,Graphics[{Red,Line[CoordinateSearchPeaks]}]]

Steepest descent.

We now try steepest descent, which requires us to know the negative gradients of our functions. We can take these gradients directly:

{D[quadF[{x,y}],x],D[quadF[{x,y}],y]}

gradQuadF[{x_,y_}] := -{24 (-2+x)+12 (-3+y),12 (-2+x)+10 (-3+y)};

The gradient of the peaks function is so unpleasant I don’t want to write it out twice!

gradPeaks[{x_,y_}] := -{-0.174` ((2.718281828459045`)^(-1.` (((-4.52`+x))^2)-1.` (((3.29` +y))^2))) (-4.52`+x)-1.218` ((2.718281828459045`)^(-1.` (((-4.5`+x))^2)-1.` (((-3.39`+y))^2))) (-4.5`+x)-0.24` ((2.718281828459045`)^(-1.` (((-3.27`+x))^2)-1.` (((3.91` +y))^2))) (-3.27`+x)-2.14` ((2.718281828459045`)^(-1.` (((-2.17`+x))^2)-1.` (((1.83` +y))^2))) (-2.17`+x)-3.2` ((2.718281828459045`)^(-1.` (((-1.91`+x))^2)-1.` (((-1.95`+y))^2))) (-1.91`+x)+0.758` ((2.718281828459045`)^(-1.` (((-1.87`+x))^2)-1.` (((1.85` +y))^2))) (-1.87`+x)+0.666` ((2.718281828459045`)^(-1.` (((-1.68`+x))^2)-1.` (((-1.63`+y))^2))) (-1.68`+x)-0.95` ((2.718281828459045`)^(-1.` (((-1.31`+x))^2)-1.` (((-0.689`+y))^2))) (-1.31`+x)-2.88` ((2.718281828459045`)^(-1.` (((-1.28`+x))^2)-1.` (((-0.8`+y))^2))) (-1.28`+x)+0.514` ((2.718281828459045`)^(-1.` (((-1.26`+x))^2)-1.` (((-2.22`+y))^2))) (-1.26`+x)-1.868` ((2.718281828459045`)^(-1.` (((-1.22`+x))^2)-1.` (((0.49` +y))^2))) (-1.22`+x)+0.292` ((2.718281828459045`)^(-1.` (((-0.95`+x))^2)-1.` (((0.716` +y))^2))) (-0.95`+x)-1.832` ((2.718281828459045`)^(-1.` (((-0.79`+x))^2)-1.` (((-4.76`+y))^2))) (-0.79`+x)-3.648` ((2.718281828459045`)^(-1.` (((-0.35`+x))^2)-1.` (((-4.32`+y))^2))) (-0.35`+x)+0.1` x-3.598` ((2.718281828459045`)^(-1.` (((0.42` +x))^2)-1.` (((3.` +y))^2))) (0.42` +x)-1.65` ((2.718281828459045`)^(-1.` (((1.34` +x))^2)-1.` (((3.742` +y))^2))) (1.34` +x)-3.36` ((2.718281828459045`)^(-1.` (((1.38` +x))^2)-1.` (((4.91` +y))^2))) (1.38` +x)-2.874` ((2.718281828459045`)^(-1.` (((1.45` +x))^2)-1.` (((2.88` +y))^2))) (1.45` +x)+0.68` ((2.718281828459045`)^(-1.` (((1.54` +x))^2)-1.` (((-1.67`+y))^2))) (1.54` +x)-3.14` ((2.718281828459045`)^(-1.` (((1.58` +x))^2)-1.` (((3.72` +y))^2))) (1.58` +x)-0.42` ((2.718281828459045`)^(-1.` (((2.32` +x))^2)-1.` (((-1.69`+y))^2))) (2.32` +x)-3.16` ((2.718281828459045`)^(-1.` (((2.36` +x))^2)-1.` (((-4.56`+y))^2))) (2.36` +x)-0.284` ((2.718281828459045`)^(-1.` (((3.13` +x))^2)-1.` (((-3.85`+y))^2))) (3.13` +x)-0.572` ((2.718281828459045`)^(-1.` (((3.44` +x))^2)-1.` (((0.8` +y))^2))) (3.44` +x)-0.054` ((2.718281828459045`)^(-1.` (((3.83` +x))^2)-1.` (((-0.83`+y))^2))) (3.83` +x)-2.42` ((2.718281828459045`)^(-1.` (((3.96` +x))^2)-1.` (((-1.54`+y))^2))) (3.96` +x)+0.1254` ((2.718281828459045`)^(-1.` (((3.97` +x))^2)-1.` (((2.99` +y))^2))) (3.97` +x)-2.94` ((2.718281828459045`)^(-1.` (((4.38` +x))^2)-1.` (((0.04` +y))^2))) (4.38` +x)-3.98` ((2.718281828459045`)^(-1.` (((4.67` +x))^2)-1.` (((1.19` +y))^2))) (4.67` +x)-3.62` ((2.718281828459045`)^(-1.` (((4.98` +x))^2)-1.` (((3.2` +y))^2))) (4.98` +x),-1.832` ((2.718281828459045`)^(-1.` (((-0.79`+x))^2)-1.` (((-4.76`+y))^2))) (-4.76`+y)-3.16` ((2.718281828459045`)^(-1.` (((2.36` +x))^2)-1.` (((-4.56`+y))^2))) (-4.56`+y)-3.648` ((2.718281828459045`)^(-1.` (((-0.35`+x))^2)-1.` (((-4.32`+y))^2))) (-4.32`+y)-0.284` ((2.718281828459045`)^(-1.` (((3.13` +x))^2)-1.` (((-3.85`+y))^2))) (-3.85`+y)-1.218` ((2.718281828459045`)^(-1.` (((-4.5`+x))^2)-1.` (((-3.39`+y))^2))) (-3.39`+y)+0.514` ((2.718281828459045`)^(-1.` (((-1.26`+x))^2)-1.` (((-2.22`+y))^2))) (-2.22`+y)-3.2` ((2.718281828459045`)^(-1.` (((-1.91`+x))^2)-1.` (((-1.95`+y))^2))) (-1.95`+y)-0.42` ((2.718281828459045`)^(-1.` (((2.32` +x))^2)-1.` (((-1.69`+y))^2))) (-1.69`+y)+0.68` ((2.718281828459045`)^(-1.` (((1.54` +x))^2)-1.` (((-1.67`+y))^2))) (-1.67`+y)+0.666` ((2.718281828459045`)^(-1.` (((-1.68`+x))^2)-1.` (((-1.63`+y))^2))) (-1.63`+y)-2.42` ((2.718281828459045`)^(-1.` (((3.96` +x))^2)-1.` (((-1.54`+y))^2))) (-1.54`+y)-0.054` ((2.718281828459045`)^(-1.` (((3.83` +x))^2)-1.` (((-0.83`+y))^2))) (-0.83`+y)-2.88` ((2.718281828459045`)^(-1.` (((-1.28`+x))^2)-1.` (((-0.8`+y))^2))) (-0.8`+y)-0.95` ((2.718281828459045`)^(-1.` (((-1.31`+x))^2)-1.` (((-0.689`+y))^2))) (-0.689`+y)+0.1` y-2.94` ((2.718281828459045`)^(-1.` (((4.38` +x))^2)-1.` (((0.04` +y))^2))) (0.04` +y)-1.868` ((2.718281828459045`)^(-1.` (((-1.22`+x))^2)-1.` (((0.49` +y))^2))) (0.49` +y)+0.292` ((2.718281828459045`)^(-1.` (((-0.95`+x))^2)-1.` (((0.716` +y))^2))) (0.716` +y)-0.572` ((2.718281828459045`)^(-1.` (((3.44` +x))^2)-1.` (((0.8` +y))^2))) (0.8` +y)-3.98` ((2.718281828459045`)^(-1.` (((4.67` +x))^2)-1.` (((1.19` +y))^2))) (1.19` +y)-2.14` ((2.718281828459045`)^(-1.` (((-2.17`+x))^2)-1.` (((1.83` +y))^2))) (1.83` +y)+0.758` ((2.718281828459045`)^(-1.` (((-1.87`+x))^2)-1.` (((1.85` +y))^2))) (1.85` +y)-2.874` ((2.718281828459045`)^(-1.` (((1.45` +x))^2)-1.` (((2.88` +y))^2))) (2.88` +y)+0.1254` ((2.718281828459045`)^(-1.` (((3.97` +x))^2)-1.` (((2.99` +y))^2))) (2.99` +y)-3.598` ((2.718281828459045`)^(-1.` (((0.42` +x))^2)-1.` (((3.` +y))^2))) (3.` +y)-3.62` ((2.718281828459045`)^(-1.` (((4.98` +x))^2)-1.` (((3.2` +y))^2))) (3.2` +y)-0.174` ((2.718281828459045`)^(-1.` (((-4.52`+x))^2)-1.` (((3.29` +y))^2))) (3.29` +y)-3.14` ((2.718281828459045`)^(-1.` (((1.58` +x))^2)-1.` (((3.72` +y))^2))) (3.72` +y)-1.65` ((2.718281828459045`)^(-1.` (((1.34` +x))^2)-1.` (((3.742` +y))^2))) (3.742` +y)-0.24` ((2.718281828459045`)^(-1.` (((-3.27`+x))^2)-1.` (((3.91` +y))^2))) (3.91` +y)-3.36` ((2.718281828459045`)^(-1.` (((1.38` +x))^2)-1.` (((4.91` +y))^2))) (4.91` +y)};

Now that we’ve got the gradient functions, we can plot them

Show[quadPlot,StreamPlot[gradQuadF[{x,y}],{x,-20,20},{y,-20,20}]]

Show[peakPlot,StreamPlot[gradPeaks[{x,y}],{x,-10,10},{y,-10,10}]]

These plots are pretty convincing!

We’re now going to run the steepest descent method:

SteepestQuad = Last /@ DirectionMethod[quadF,gradQuadF,{10,-17},30];

Show[quadPlot,Graphics[{Red,Line[SteepestQuad]}]]

You can see that it’s surprisingly not much better! You can also notice the interesting fact that the successive search directions are always perpendicular to each other.

SteepestPeaks = Last /@ DirectionMethod[peaks,gradPeaks,{-4.2,-5.3},15];

Show[peakPlot,Graphics[{Red,Line[SteepestPeaks]}]]

For the peaks function, something quite similar happens.

Conjugate gradient.

So conjugate gradient is a somewhat different beast; for one thing, we need to know previous gradient values in order to compute the next conjugate direction. So it won’t quite fit in our previous DirectionMethod framework. But it’s not much longer-- we just need to add another parameter which records the current conjugate direction h and gradient g at the previous x value:

CGStep[{f_,GradFunction_,x_,g_,h_}] :=Module[{v,t,γ,ξ},
t=FindMinimum[{f[x + t h],t > 0},{t}][[2,1,2]];
ξ=GradFunction[x + t h];
γ=-(((-ξ+g).ξ)/(g.g));
{f,GradFunction,x+t h,ξ,ξ + γ h}
];

CGMethod[f_,Dir_,x0_,n_] := NestList[CGStep,{f,Dir,x0,Dir[x0],Dir[x0]},n];

CGQuad = #[[3]]& /@ CGMethod[quadF,gradQuadF,{10,-17},2]

Show[quadPlot,Graphics[{Red,Line[CGQuad]}]]

CGPeaks = #[[3]]& /@ CGMethod[peaks,gradPeaks,{-4.2,-5.3},8]

Show[peakPlot,Graphics[{Red,Line[CGPeaks]}],Graphics[{Green,Line[SteepestPeaks]}]]

Show[ContourPlot[peaks[{x,y}],{x,-1.5,-1},{y,1,1.5},Contours->35],Graphics[{Red,Line[CGPeaks]}],Graphics[{Green,Line[SteepestPeaks]}],PlotRange->{{-1.5,-1},{1,1.5}}]