The Nelder-Mead Simplex Algorithm for Numerical Minimization.

In this lecture, we'll explore the Nelder-Mead method for function optimization. This is a relatively poorly understood, but incredibly useful and commonly used method for minimizing functions of many variables when you don't have derivatives available. 

Lagarias has shown that the Nelder-Mead algorithm converges for strictly convex functions on R^2 (the paper is posted on the course webpage), but the proof is too intricate for us to do in class (I think). So instead we'll follow the treatment in Matthews/Fink "Numerical Methods Using MATLAB", which is also posted on the course page.

We saw in the class discussion that there are four basic move types:

Reflection[B_,G_,W_] := W + 2((B + G)/2  - W);
Expansion[B_,G_,W_] := W + 3((B + G)/2 - W);
Contraction[B_,G_,W_,R_,f_] := Module[{M,C1,C2},
M = (B+G)/2;
C1 = (M + W)/2;
C2 = (M + R)/2;
If [f[C1] < f[C2],C1,C2]
];
Shrink[B_,G_,W_,f_] := If[f[(B+W)/2] < f[(B+G)/2],{{B,(B+W)/2,(B+G)/2},f},{{B,(B+G)/2,(B+W)/2},f}];

We can now write a short piece of code which encapsulates the Nelder-Mead logic from the class lecture.

NMStep[{{BB_,GG_,WW_},f_}] := Module[{B,G,W,R,E,C},

{B,G,W} = SortBy[{BB,GG,WW},f];

R = Reflection[B,G,W];
	
If [f[R] < f[G],

	If[f[B] < f[R],{{B,G,R},f},
	
	E = Expansion[B,G,W];

	If[f[E] < f[B], {{B,G,E},f}, {{B,G,R},f}] ]

	,
	
	If [f[R] < f[W],{{B,G,R},f},

	C = Contraction[B,G,W,R,f];

	If [f[C] < f[W],{{B,G,C},f},Shrink[B,G,W,f]]

	]

]

];

We now use this to put together an interactive demo for Nelder-Mead in just a few lines of code! We could make this demo more efficient by computing the entire set of Nelder-Mead steps in advance using NestList, but that would be more complicated (and anyway, this is fast enough for our purposes).

Clear[f];

NelderMeadDemo[f_] := Module[{B,G,W,x,y,pt1,pt2,pt3},
fplot = ContourPlot[f[{x,y}],{x,-5,5},{y,-5,5},Contours->40,PlotRange->Full];
{B,G,W} = SortBy[{{1.0,2.0},{2.0,4.0},{3.0,2.0}},f];
Manipulate[
If[i == 0,
Show[fplot,Graphics[{FaceForm[Opacity[0.2]],EdgeForm[Thick],Polygon[{pt1,pt2,pt3}]}],ImageSize->600],
Show[fplot,Graphics[{FaceForm[Opacity[0.2]],EdgeForm[Thick],Polygon[Nest[NMStep,{{pt1,pt2,pt3},f},i][[1]]]}],ImageSize->600]
],
{{pt1,B},Locator},{{pt2,G},Locator},{{pt3,W},Locator},{{i,0,"Number of Nelder-Mead Steps"},0,20,1}]
]

We now define a couple of test functions, including a version of the "peaks" function from our demo on direction search methods.

f[{x_,y_}] := x^2 + x*Sin[y] + 0.5* y^2;

PeakLocs = {{1.2645164185495439,2.2209592594231644},{3.270097950964729,-3.9103464721409935},{-1.5418494655660115,1.67768310221404},{1.2801691690766894,0.8052153454360766},{-1.5862482581167558,-3.7255065039308306},{-1.4525381551315597,-2.889011754305841},{0.9500722910632486,-0.7161119376122151},{1.312202342159729,0.6899999139503841},{4.503567692804019,3.395568485639089},{-2.369718842541616,4.562190074588269},{1.6845592922341481,1.6356335849355972},{-3.4468198549466056,-0.800970903025604},{1.22112277017537,-0.49751814087780133},{-3.970165263007001,-2.99289767452934},{-0.4208656595723781,-3.001849169975772},{-1.3856312348971755,-4.915760721595204},{-3.968840948897043,1.5461551708825247},{-4.384214376606435,-0.047803586643848206},{-3.134811414251807,3.8566937251845204},{1.9170927826948194,1.9515643677788574},{0.35504210470119446,4.321129208956641},{-2.326302192066418,1.696388350640209},{-3.835147627998543,0.8335337362073698},{-4.677050568662207,-1.1921397281305666},{1.8701748149086974,-1.8588104922747934},{4.529280893642241,-3.29069439562174},{-1.3413052778439631,-3.7426290421130144},{-4.988037916339003,-3.205943557916985},{2.1751817275642935,-1.831364822604339},{0.7939885412091634,4.760878143088036}};

PeakHeights ={-0.25777991514805754,0.12168890305875868,-0.3401322490537457,1.448763518105869,1.571140898920162,1.4374031291880351,-0.14601876542537306,0.4756815200321389,0.6099270368096716,1.5823576233024896,-0.33339080193342596,0.28696927385167204,0.9348230176039847,-0.06278226492318018,1.7998741291598934,1.6800123227287589,1.2122844377642101,1.4741373935483724,0.14204449977632683,1.6000531442041135,1.824457699684281,0.21097984758539168,0.027757089970823934,1.999020600768585,-0.3792507004414136,0.08793575013506594,0.825554512431145,1.8165403432053129,1.07144672788273,0.9166011753859027};

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

NelderMeadDemo[f]

Next time, we'll have a discussion of the Nelder-Mead algorithm in high dimensions, with a corresponding high-dimensional demo!