One-Dimensional Minimization

We are interested in finding the minimum of a function without using derivative information.

Golden Section Minimization

We have already discussed the golden section search algorithm (Fibonacci search is just too much of a pain to code!). 

We recall that it works in the following way:

f[x_] := Sin[x];

Golden[{a_,b_}] := Module[{r},

r = 1.0/(0.5 + Sqrt[5]); (* Force machine arithmetic *)
x = a + r*(b-a);
y = a + r^2*(b-a);

If [f[x] < f[y],{N[y,15],N[b,15]},{N[a,15],N[x,15]}]

];

We can now generate a list of Golden section search points using NestList.

GoldenPoints = NestList[Golden,{4.0,8.0},100]

DrawInterval[{a_,b_}] := Module[{r,xx,yy},

r = 1.0/(0.5 + Sqrt[5]); (* Force machine arithmetic *)
xx = a + r*(b-a);
yy = a + r^2*(b-a);

{Plot[f[x],{x,a,b},Filling->Axis,FillingStyle->Automatic,AxesOrigin->{0,0}],
ListPlot[{{a,f[a]},{b,f[b]},{a,0},{b,0}},PlotMarkers-> {Automatic,Small}],
Plot[f[x],{x,xx,yy},Filling->Axis,FillingStyle->Automatic,AxesOrigin->{0,0}],
ListPlot[{{xx,f[xx]},{xx,0},{yy,f[yy]},{yy,0}},PlotMarkers-> {Automatic,Small}]
}

]

Animate[Show[Plot[f[x],{x,4.0,5.0},Filling->Axis,AxesOrigin->{4,0}],DrawInterval[GoldenPoints[[i]]]],{i,1,20,1}]

Of course, we know that the minimum is at 3 pi/2, so we can compute the number of correct digits as usual.

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];

GoldenDigits = Table[Digits[(GoldenPoints[[i,2]] + GoldenPoints[[i,1]])/2.0,3*Pi/2],{i,1,Length[GoldenPoints]}]

ListLinePlot[{GoldenDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

Now it's interesting to note that Golden will continue to subdivide the interval until the size of the interval reaches machine epsilon. (It just makes the wrong choices about how to do it.) So if we plot the width of the interval as an approximation to zero, we get

GoldenWidth= Table[Digits[(GoldenPoints[[i,2]]- GoldenPoints[[i,1]]),0],{i,1,Length[GoldenPoints]}]

ListLinePlot[{GoldenWidth},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

So don't be fooled by a very small interval coming back from a search of this type. The actual minimum is probably not in the interval in the first place!

Quadratic Interpolation.

In quadratic interpolation, we fit a quadratic polynomial to our function at three points and minimize the polynomial. Deriving the algorithm is the sort of thing that's most easily done with Mathematica:

Clear[f]; Clear[x];

IP = InterpolatingPolynomial[{{a,f[a]},{b,f[b]},{c,f[c]}},x]

NextPoint = x /. Simplify[Solve[D[IP,x] == 0,x]]

MinTest = Simplify[D[IP,x,x]]

We now have a function which gives us our next guess in terms of the previous three points and we have a test to see whether this is expected to be a minimum or a maximum. 

We can now construct an algorithm. Observe that you can easily use Mathematica to write Mathematica code by translating your output into // InputForm.

NextPoint // InputForm

MinTest // InputForm

Here's the basic idea of the algorithm. Suppose we're given three points in which to search for the min.

0. If the three points are close enough together, we don't try to evaluate anything else.  If not...

1. Compute the minimum of the quadratic interpolating polynomial and the second derivative of the interpolating polynomial at that point.

2. If the interpolating polynomial has a min at the point in question, 

	if it is close enough to our triple, we take it
	if not, we step in that direction.
	
3. If the interpolating polynomial has a max at the point in question, we just try to step ``downhill''.

Tol = 10^(-6);
MaxStep = 10^(-2);

QuadraticInterpolation[{aa_,bb_,cc_}] := Module[{qMin,deriv,distFromABC,a,b,c},

{a,b,c} = Sort[{aa,bb,cc}]; 

If[ c -a < Tol,Return[{a,b,c}],

qMin =((-b^2+c^2)*f[a]+(a^2-c^2)*f[b]+(-a^2+b^2)*f[c])/(2*((-b+c)*f[a]+(a-c)*f[b]+(-a+b)*f[c]));
	
deriv = (2*((-f[a]+f[b])/(a-b)+(f[b]-f[c])/(b-c)))/(-a+c);
	
If[ deriv > 0, 

(* The new point is a min *)

Which[ (* This evaluates to the first True conditional *)

qMin < a,
If [a - qMin < MaxStep, {qMin,a,b},{a - MaxStep,a,b}],

qMin < b,
{a,qMin,b},

qMin < c,
{b,qMin,c},

qMin >= c,
If [qMin - c < MaxStep, {b,c,qMin},{b,c,c+MaxStep}]

],

(* The new point is a max *)

If [f[a] < f[c], {a - MaxStep,a,b}, {b,c,c+MaxStep}]

] 

]

];

f[x_] := Sin[x];
QuadraticInterpolationPoints = NestList[QuadraticInterpolation,{4.0,5.0,6.0},100]

QIDigits = Table[Digits[(QuadraticInterpolationPoints[[i,2]] + QuadraticInterpolationPoints[[i,1]])/2.0,3*Pi/2],{i,1,Length[QuadraticInterpolationPoints]}]

ListLinePlot[{GoldenDigits,QIDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

What's happening here? Essentially, when the points get too close together, the quality of our quadratic interpolation goes south due to numerical error! 

If we ask for a smaller Tol, we can stablize the situation (at a lower accuracy). 

The essential problem here is that we gave up the certainty of having the minimum bracketed within the interval [a,c] in order to take large quadratic interpolation steps.

Brent's Method for Minimization

The idea of Brent's method is that we should alternate between quadratic minimization and golden section steps, using the fast convergence of quadratic minimization when it's working and the guaranteed convergence of golden section minimization when it's not. 

Here, then, is a version of Brent's "localmin", translated from the original Algol 60 into Mathematica. The basic outline of the algorithm is (recall) that 

a,b : an interval that contains the min.
x    : the location of the smallest f value found yet
w   : the location of the next smallest f value found yet
v    : the location of the third smallest value of f found yet
u    : the last point where f was evaluated (regardless of the
        value found for f at that point)
        
At each stage, we fit a parabola through values of f at x, v, and w. If the parabolic fit obeys certain conditions (it's a min, the step size is decreasing fast enough, it's in the search interval) we take it. Otherwise, we do a golden section step. We stop when the search interval is small enough. 

We deal with roundoff error by supplying two parameters: eps and t. The error tolerance on a given step should be

tol = eps*x + t.

For the reasons discussed above, eps should be something like Sqrt[machine epsilon] and t should be about the same.

BrentMin[aa_,bb_,eps_,t_,f_] := Module[{a,b,c,d,e,m,p,q,r,tol,u,v,w,fu,fv,fw,fx,x,IntList,StepCount},

(* These globals will be reset every time we run BrentMin. *)

QuadraticSteps = 0;
GoldenSteps = 0;

(* We start by initializing the variables *)

a = 1.0*aa;
b = 1.0*bb;
c = (3.0 - Sqrt[5.0])/2.0;

v = a + c * (b-a);
w = a + c * (b - a);
x = a + c * (b-a);

(* We will later compute a quadratic fit step p/q. e is the second to last value of this step. Since we havent stepped yet, we initialize e now. *)

e = 0;                               

fx = f[x];
fv = f[v];
fw = f[w];

m = (0.5)*(b+a); (* the midpoint of [a,b] *)
tol = eps*Abs[x] + t;
Print[tol];

IntList = {{a,x,b}}; (* list of working intervals *)
StepCount = 0;

Print["Initial interval:",IntList];
(*The stopping criterion is whether x is too close to *)
(* midpoint*)

While[Abs[x - m] > (2*tol - 0.5*(b-a)) && StepCount < 40,

	Print["{a,x,b} = ",{a,x,b}// InputForm ];
	IntList = Append[IntList,{a,x,b}];
	StepCount++;

	m = (0.5)*(b+a);
	tol = eps*Abs[x] + t; 
	p = 0; q = 0; r = 0;

	If[ Abs[e] > tol ,

(* The parabolic fits are still good enough for us to *)
 (*  try again. We now compute a new fit. This is the *)
   (* code for the minimum of the quadratic polynomial *)
(* through x, v, and w, which will be at x + p/q. *)

		r = (x - w)(fx - fv); 
		q = (x - v)(fx - fw);
		p = (x - v)  q  - (x - w) r; 
		q = 2.0 (q - r);

		If [ q > 0, p = -p, q = -q];

		r = e; e = d;

	];

	(* We want to check whether 

	p/q <(1/2) r 

	without actually dividing by q (in case q is zero). 
   So instead we check p < (1/2) qr. 

   We also need to check whether 

	a < x + p/q < b

	again, without actually dividing by q. We get *)

	(*Print["p = ", p,"q = ",q," Abs[0.5*q*r] = ",Abs[0.5 * q * r], " p = ",p," q(a-x) = ",q(a-x),"q (b-x) = ", q(b-x),"Res = ", Abs[p] < Abs[0.5 * q * r] ∧ p > q(a-x) ∧ p <q(b-x)]; *)

	If [ Abs[p] < Abs[0.5 * q * r] ∧ p > q(a-x) ∧ p <q(b-x),

	(* Now we're going to do a parabolic interpolation step. *)
	(* we won't get here if q = 0. *)

	d = p/q;
	u = x + d;

	Print["\t Quadratic step to ",u];
	QuadraticSteps++;
       (* We now make sure that u is far enough from a and b*)

	      If [ u - a < 2 tol || b - u < 2 tol, d = If [x < m,tol,-tol]];

          ,

	     (* Do a golden section step *)

	       e = If[ x < m, b, a] - x; d = c e;
	       Print["\t Golden step to ",x+d];
	       GoldenSteps++;

	  ];

	(* We now want to compute u and update fu. But we need to make 
      sure that u is not too close to x *)

	u = x + If[ Abs[d] ≥ tol, d, If[d > 0,tol,-tol]];
	fu = f[u];

	(*Print["x = ",x," fx = ",fx," u = ",u," fu = ",fu," w = ",w," fw = ",fw,"fu - fx = ",fu - fx, "fu <= fx = ", fu ≤ fx,"u - x= ",u-x,"u < x =", u < x]; *)

	(* Now we need to update the search interval and the list of 
	  best, second best, and third best locations x, w, v *)

	If [ fu ≤ fx,
	(* We use the pair [u,x] to update interval to [a,x] or [x,b] *)

	(*Print["fu <= fx"];*)
	If[u<x,b=x,a=x];

	(* Now we update the rest of the interval *)
	v = w; fv = fw; w = x; fw = fx; x = u; fx = fu;
	,
	(* We use the pair [u,x] to update interval to [a,u] or [u,b] *)
	If[u<x,a=u,b=u];
	
	(* fu is not < fx, so u does not replace x. Does it replace w? *)
	If[fu≤fw || w == x,
	(*yes, so change v, too*)
	v=w; fv = fw; w = u; fw = fu;
	,
	(*no, so check whether it replaces v *)
	If[fu ≤ fv || v == x || v ==w ,v = u; fv = fu;]

       ](* end of if fu < fw *)

      ] (* end of if fu < fx *)

     (* At this point, we've reset the interval, so we store the data. *)

];(* end of main while loop *)

IntList (* return the list of trial intervals *)

];

BrentIntervals = BrentMin[4,8,10.0^(-7),10.0^(-7),Sin ]

BrentDigits = Table[Digits[BrentIntervals[[i,2]],3.0*Pi/2.0],{i,1,Length[BrentIntervals]}];
BrentDigits = PadRight[BrentDigits,100,Max[BrentDigits]]

ListLinePlot[{GoldenDigits,QIDigits,BrentDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

We can see from this the tremendous advantage of Brent's method! Even though we only deserve 8 digits of accuracy, we actually get 10. How is this possible? Well, recall that we're converging through parabolic approximation to the min as long as possible, and that method converges superlinearly (Brent gives conditions for this to converge with exponent 1.3427..., and Sin certainly meets these criteria). 

So we see that the number of correct digits doubles on the second to last step (from 5 to 10), as we take a really good quadratic fit to our function.

More difficult functions

Of course, Sin[x] is not a very difficult function to minimize. Bad examples include functions that are flatter than quadratic at the minimum, such as x^4.

f[x_] := x^4;
GoldenPoints = NestList[Golden,{-1,3},100]
GoldenDigits = Table[Digits[(GoldenPoints[[i,2]] + GoldenPoints[[i,1]])/2.0,0],{i,1,Length[GoldenPoints]}];
BrentIntervals = BrentMin[-1,3,10.0^(-7),10.0^(-7),f ] 
BrentDigits = Table[Digits[BrentIntervals[[i,2]],0],{i,1,Length[BrentIntervals]}];
BrentDigits = PadRight[BrentDigits,100,Max[BrentDigits]];

Length[BrentIntervals]

ListLinePlot[{GoldenDigits,BrentDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

We see that Brent's method terminates after 15 function evaluations. Golden section search continues to chug along, though! By blind luck, it manages to make the right decisions on the min and eventually catches up after 60 evaluations or do.

Another weird example comes when you consider functions that aren't very smooth:

f[x_] := Piecewise[{{-x,x<-1},{x^2 - Cos[x*Pi/2],x<1},{x,x>1}}];
Plot[f[x],{x,-2,2}]

GoldenPoints = NestList[Golden,{-1,3},100]
GoldenDigits = Table[Digits[(GoldenPoints[[i,2]] + GoldenPoints[[i,1]])/2.0,0],{i,1,Length[GoldenPoints]}];
BrentIntervals = BrentMin[-1,3,10.0^(-7),10.0^(-7),f ] 
BrentDigits = Table[Digits[BrentIntervals[[i,2]],0],{i,1,Length[BrentIntervals]}];
BrentDigits = PadRight[BrentDigits,100,Max[BrentDigits]];

ListLinePlot[{GoldenDigits,BrentDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

Functions which aren't smooth at the minimum can cause both golden section search and Brent's method to struggle a bit:

f[x_] := Piecewise[{{x^4,x<-1},{Abs[x] ,x<1},{x^3,x>1}}];
Plot[f[x],{x,-2,2}]

GoldenPoints = NestList[Golden,{-1,3},100]
GoldenDigits = Table[Digits[(GoldenPoints[[i,2]] + GoldenPoints[[i,1]])/2.0,0],{i,1,Length[GoldenPoints]}];
BrentIntervals = BrentMin[-1,3,10.0^(-7),10.0^(-7),f ] 
BrentDigits = Table[Digits[BrentIntervals[[i,2]],0],{i,1,Length[BrentIntervals]}];
BrentDigits = PadRight[BrentDigits,100,Max[BrentDigits]];

ListLinePlot[{GoldenDigits,BrentDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]