Minimization with Derivatives

We are still interested in finding the minimum value of a function of one variable, but this time we assume that we have one derivative of the function as well. 

We are going to use a variant of Brent's method (originally proposed by Brent, and spelled out in detail in the excerpt from Numerical Recipes on the class webpage).

dBrentMin[ax_,bx_,cx_,f_,df_,tol_] := Module[{a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm,iter,ok1,ok2,IntList},

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

InterpolationSteps = 0;
BisectionSteps = 0;
dBrentFEvals = 0;
dBrentdFEvals = 0;
iter = 0;

(* We start by initializing the variables *)

e = 0.0;

a = If[ax < cx, ax, cx];
b = If[ax > cx, ax, cx];

x = bx; w = bx; v = bx;
fw = f[x]; fv = f[x]; fx = f[x]; dBrentFEvals++;
dw = df[x]; dv = df[x]; dx = df[x];dBrentdFEvals++;

IntList = {{a,x,b}}; (* list of working intervals *)
Print["Initial interval:",IntList];

For[iter=1,iter<100,iter++, (* This will loop, counting Iter, as long as Iter < 100 *)

	xm = 0.5*(a+b);
	tol1 = tol*Abs[x]+10^(-10);
	tol2 = 2.0*tol1;

	If [Abs[x-xm] ≤ (tol2 - 0.5*(b-a)), (*the termination criterion*)
	Print["Returning",{x,fx}," (interval < tol)"]; Return[{x,f[x]}], a = a];  (* the a=a just does nothing *)

	If[Abs[e] > tol1,  (* e is the movement on the last step *)

	d1 = If[dw ≠ dx, (w-x)*dx/(dx - dw), 2.0*(b-a)];
	d2 = If[dv ≠ dx, (v-x)*dx/(dx - dv), 2.0 *(b-a)];

	a=a; (* We guess using secant method on df at x,w and x,v *)

	u1 = x + d1;  (* These are the trial positions given by secant *)
	u2 = x + d2; 	

	ok1 = TrueQ[(a - u1)*(u1 -b) > 0.0 && dx*d1 ≤ 0.0];
	ok2 = TrueQ[(a - u2)*(u2 - b)> 0.0 && dx*d2 ≤ 0.0];

	olde=e;
	e = d;

	a=a; (* Check that each u is in the interval and that the derivative at x is pointing downhill towards u. *)

	If[ ok1 || ok2, (* If one of these is ok, take it (or the smallest if both are ok *)
	
	d = If [ ok1 && ok2, If [Abs[d1] < Abs[d2], d1, d2], If[ok1, d1, d2]];
	If[Abs[d] ≤ Abs[0.5*olde], (* check for sufficient progress *)
	u=x+d; 
	d = If[u-a < tol2 || b-u < tol2, (* check u too close to endpt *)
	If[xm-x>0.0,tol1,-tol1],d];
	Print["Interpolation step"];
	,
	Print["Bisection step (interpolation not converging) dx = ",dx];
	e = If[dx >= 0.0, a-x,b-x]; (* bisect (d too big) *)
	d = 0.5*e; 
	];
	,
	  Print["Bisection step (interpolation points outside [a,b]) dx = ",dx];
	  e = If[dx >= 0.0,a-x,b-x]; (* bisect (u1, u2 not ok) *)
	  d = 0.5*e;
       ];
,
Print["Bisection step (interpolation steps are too small) dx =",dx];
e = If[dx >= 0.0,a-x,b-x]; (* bisect (e ≤ tol1) *)
d = 0.5*e;	
];

If[Abs[d] ≥ tol1, (* make sure d is large enough *)
u = x + d;
fu = f[u]; dBrentFEvals++;
, 
u = x + If[ d >0,tol1,-tol1]; (* if not, replace d by tol1 *)
fu = f[u]; dBrentFEvals++;
If [fu > fx || (a - u)*(u -b) <0.0, Print["Returning " ,{x,fx}," (min step went up or out)"]; Return[{x,fx}],u=u]; (* if the min step goes up or out of the interval, quit *)
];

du = df[u]; dBrentdFEvals++;

a=a; (* Now we need to update the interval and x, w, v etc. *)

If [fu ≤ fx, (* the interval will now be [a,x] or [x,b] *)
	
	If [ u ≥ x, a = x, b =x];
	
	{v,fv,dv} = {w,fw,dw}; (* u replaces x *)
	{w,fw,dw} = {x,fx,dx};
	{x,fx,dx} = {u,fu,du};

, (* the interval will now be [a,u] or [u,b] *)

	If [u < x, a = u, b = u];

	If [fu ≤ fw || w ==x, (* u replaces w *)

	{v,fv,dv} = {w,fw,dw};
	{w,fw,dw} = {u,fu,du};

	, 

	If [fu < fv || v ==x, (* u replaces v *)
		
		{v,fv,dv} = {u,fu,du},a = a];

	];

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

];

Print["Too many iterations",iter];
IntList

];

We now try out this method on the standard Sin[x] problem on [4,8]. We'll choose 5 as the middle point of the bracketing triple. We'll let 10^(-7) be the tolerance. Note that we're keeping count of the function evaluations and the derivative evaluations as we go.

dBrentMin[4.0,5.52,8.0,Sin,Cos,10^(-7)]
dBrentFEvals
dBrentdFEvals

Comparing dBrentMin to BrentMin.

We now want to compare these results to the standard Brent's method. We modify the code slightly to count function evaluations in BrentFEvals and to return a final answer in the same form as dBrentMin.

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;
BrentFEvals = 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);

 a= 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]; BrentFEvals++;
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 ,

a=a;(* The parabolic fits are still good enough for us to *)
a=a;(*  try again. We now compute a new fit. This is the *)
a=a;(* code for the minimum of the quadratic polynomial *)
a=a;(* 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++;
       a =a; (* 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]; BrentFEvals++;

	(*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 *)

x = RotateRight[IntList ][[1,2]];(* return the list of trial intervals *)
{x,f[x]}
];

BrentMin[4.0,8.0,10^(-7),10^(-7),Sin]
BrentFEvals

We see here an interesting piece of information! While the dBrent method converged in fewer steps, it did not actually use fewer function evaluations: each method used a total of 11 function (or derivative) evaluations. Since a function evaluation is about as costly as a derivative evaluation in general, this shows that dBrent doesn't always improve performance (at least on relatively easy problems like this one).

More example problems for BrentMin and dBrentMin

We now try some more examples for BrentMin and dBrentMin. Of course, dBrentMin is not going to be applicable on the problems which are not differentiable. But we can still try functions which have no second derivatives. We first try a very flat function.

f[x_] := (x+1)^4  + 12( x+1)^2 + 1;
df[x_] := D[f[u],u] /. {u -> x};
Plot[f[x],{x,-3,1}]

dBrentMin[-3.0,-2.0,12.0,f,df,10^(-7)]
dBrentFEvals + dBrentdFEvals

BrentMin[-3.0,12.0,10^(-7),10^(-7),f]
BrentFEvals

Observe that there's no time advantage to dBrent (in fact, it's slower, with 21 function evals vs 13 for Brent), but it does give you a better answer. 

We now want to try some functions which are differentiable, but not twice differentiable. It's not so hard to get that by fitting polynomials and trig functions together. Notice that this kind of stuff is more troubling for the more delicate dBrent solver than for the traditional Brent's method!

f[x_] := Piecewise[{{-Cos[x+1],x < -1},{3*(x+1)^2 -1, x≥ -1}}];
df[x_] := Piecewise[{{Sin[x+1], x < -1}, {6*(x+1), x ≥ -1}}];
Plot[f[x],{x,-3,1}]

dBrentMin[-3.0,-2,5,f,df,10^(-7)]
dBrentFEvals + dBrentdFEvals

BrentMin[-3.0,5,10^(-7),10^(-7),f]
BrentFEvals

What if we tried this with numerical derivatives?

Suppose that instead of using Cos[x] for the derivative of Sin, we picked an h, say our best value of 10^(-5) and evaluated a numerical approximation to the derivative instead using the central difference formula

f'(x) = ((f(x+h) - f(x-h))/(2h))

We get

f[x_] := Sin[x];
df[x_] := Module[{h},
h = 10.0^(-5.0);
(Sin[x+h] - Sin[x - h])/(2*h)
];

dBrentMin[4.0,5.0,8.0,f,df,10^(-7)]
dBrentFEvals+2*dBrentdFEvals

BrentMin[4.0,8.0,10^(-7),10^(-7),f]
BrentFEvals

Even though this actually does converge (roundoff error didn't get us yet!), you can't beat the dealer with this sort of thing: dBrent already uses as many function + derivative evaluations as Brent uses function evaluations, and numerically approximating the derivative means evaluating 2 (or for Richardson extrapolation, 4 or more!) function points for each derivative approximation. 

Much better to just use Brent and be done with it.