Fletcher's Algorithm for Inexact Line Minimization.

We now want to demonstrate the use of the modified Goldstein conditions in inexact line search. We're following the description of the algorithm in Fletcher "Practical Methods of Optimization" (2nd edition, pages 34 and 35).

The basic idea is first to bracket the approximate minimizer by starting at alpha = 0 and jumping out to the right by progressively larger intervals. 

We will specify in advance a value L for f small enough that we're willing to accept any alpha value with f(alpha) < L. Eventually, we'll either get to such an f(alpha) or we'll generate a bracketing interval. 

Fletcher proves (Lemma 2.6.1) that such an interval will contain a subinterval of points satisfying the modified Goldstein conditions: we then find an inexact minimizer in that interval by using a sectioning algorithm.

The core of this method will involve finding the minimum of a cubic polynomial interpolating at points a and aOld the function values fa,faOld and the derivative dfa, dfaOld. We can work out this formula with Mathematica:

IP = InterpolatingPolynomial[{{a,fa,dfa},{b,fb,dfb}},x]

Simplify[CritPts = Solve[D[IP,x] == 0, x]] // InputForm

Simplify[D[IP,x,x] /. CritPts] // InputForm

IP2 = InterpolatingPolynomial[{{a,fa,dfa},{b,fb}},x]

Simplify[CritPts = Solve[D[IP2,x] == 0, x]] // InputForm

We notice something neat about these formulae: the second derivative is in the form +/- 2 * Sqrt[stuff]/(a-b)^3. This means that the critical point which is the min (second derivative positive) is the first critical point if b > a and the second if b < a. We can then write:

CubicMin[{{a_,fa_,dfa_},{b_,fb_,dfb_}}] := Module[ {critpt1,critpt2,dd1, dd2},

If[(3*(a*(dfa+dfb)-b*(dfa+dfb)-2*fa+2*fb))< 10^(-10), (* Our cubic is actually quadratic! In this case, the formulae below will crash.*)
If[ (2*a*dfa-2*b*dfa-2*fa+2*fb)  < 10^(-10),
(* Our cubic is actually linear! Return the lower endpoint. *)
 If [fa < fb, fa,fb],
 (a^2*dfa-b^2*dfa-2*a*fa+2*a*fb)/(2*a*dfa-2*b*dfa-2*fa+2*fb)
],
If [b> a,
-(-(a^2*dfa)-a*b*dfa+2*b^2*dfa-2*a^2*dfb+a*b*dfb+b^2*dfb+3*a*fa+3*b*fa+Sqrt[(a-b)^2*(a^2*(dfa^2+dfa*dfb+dfb^2)+b^2*(dfa^2+dfa*dfb+dfb^2)-2*a*(b*(dfa^2+dfa*dfb+dfb^2)+3*(dfa+dfb)*(fa-fb))+6*b*(dfa+dfb)*(fa-fb)+9*(fa-fb)^2)]-3*a*fb-3*b*fb)/(3*(a*(dfa+dfb)-b*(dfa+dfb)-2*fa+2*fb))
,
 If[b<a,
(a^2*dfa+a*b*dfa-2*b^2*dfa+2*a^2*dfb-a*b*dfb-b^2*dfb-3*a*fa-3*b*fa+Sqrt[(a-b)^2*(a^2*(dfa^2+dfa*dfb+dfb^2)+b^2*(dfa^2+dfa*dfb+dfb^2)-2*a*(b*(dfa^2+dfa*dfb+dfb^2)+3*(dfa+dfb)*(fa-fb))+6*b*(dfa+dfb)*(fa-fb)+9*(fa-fb)^2)]+3*a*fb+3*b*fb)/(3*(a*(dfa+dfb)-b*(dfa+dfb)-2*fa+2*fb))
,
a]
]
]

];

Of course, it would be highly unwise to trust this formula without some sanity checks:

TestPoly = x^3 + 12x^2 + 14*x - 3;
dTestPoly = D[TestPoly,x]
Plot[TestPoly,{x,-10,5}]
data = {{-4,TestPoly /. x -> -4, dTestPoly /. x-> -4},{2,TestPoly /. x-> 2, dTestPoly /. x-> 2}}
N[CubicMin[data]]
N[CubicMin[RotateRight[data,1]]]

The Core of the Fletcher Minimization Algorithm

The core of the Fletcher algorithm is the routine for producing an initial bracket. This should either terminate with an acceptable point a or with an interval {a,b} (we don't know a < b) so that 

a is the best (least value of f) trial point so that f(a) < f(0) + a rho f'(0).
f'(a) satisfies (b - a) f'(a) < 0 and |f'(a)| > -sigma f'(0).
b has either f(b) > f(0)  + b rho f'(0) or f(b) > f(a) or both

Under these conditions, {a,b} is guaranteed to contain an interval of points that satisfy the modified Goldstein conditions. (Note that a violates the second modified Goldstein condition-- if it didn't, we'd simply return it here and forget about the bracketing step of the algorithm).

FletcherMin[f_,df_,L_] := Module[{rho = 0.01,sigma = 0.9, tau1 =9.0,mu,a,aOld,fa,faOld,dfa,dfaOld,fzero, dfzero, iter,targetA,targetB,alpha,alphaOld},

If[df[0] == 0,Return[{0,f[0]}]];
If[df[0] > 0,Print["df[0] > 0, Fletcher's method doesn't apply"]; Return[];];

aOld = 0;
faOld = f[aOld]; dfaOld = df[aOld]; 	
fzero = f[0]; dfzero = df[0];

mu = (L - fzero)/( rho*dfzero); 
a = Min[{mu/10.0,1.0}];(* Must choose in [0,mu] *)
fa = f[a]; dfa = df[a]; 

FEvals= 3; dFEvals=3;

For[iter=1,iter<100,iter++,

	Prinp["a = ",a];
	fa = f[a]; FEvals++;
	
	If [fa ≤ L, Return[{a,fa}]];
	
	If [fa >  fzero + a*rho*dfzero || fa > faOld, 
	(*Print["call1: a=",a," fa=",fa,"fzero =",fzero," a*rho*dfzero=",a*rho*dfzero,"faOld =",faOld];*) Return[FletcherBracket[{aOld,a},f,df,rho,sigma]]];

	dfa = df[a]; dFEvals++;

	If[Abs[dfa] ≤ -sigma * dfzero, (*Print["small dfa =",dfa,"-sigma*dfzero=",-sigma*dfzero];*) Return[{a,fa}]];
				    
          If[dfa ≥ 0, (*Print["call2"];*)Return[FletcherBracket[{a,aOld},f,df,rho,sigma]]];

	If [mu ≤ 2a - aOld,
	a = mu;
	,
	{targetA,targetB} = {2*a - aOld, Min[{mu,a + 9*(a - aOld)}]}; 
	a = CubicMin[{{aOld,faOld,dfaOld},{a,fa,dfa}}];
	
	(*Print["a = ",a,"aOld = ",aOld," TargetA = ",targetA]; *)
	"The next guess should be in {targetA,targetB}";
	a = If[a > targetB, targetB,a];  
	a = If [a < targetA, targetA, a]; 
	];

	If[a<0,Print["Warning: a < 0. a =",a," aOld=", aOld]];

	 {aOld,faOld,dfaOld} = {a,fa,dfa};

];

Print["Too many iterations in FletcherMin"];

];

The bracketing portion.

We now go ahead and work on the bracketing part of the algorithm. (Note that we can't test anything until the bracketing portion of the algorithm is defined.)

FletcherBracket[{aStart_,bStart_},f_,df_,rho_,sigma_] := Module[{a,b,fa,fb,dfa,dfb,tau2 = 1.0/10.0,tau3 = 1.0/2.0,targetA,targetB,j,alpha,fzero,dfzero,falpha,dfalpha},

Prinp["Starting FletcherBracket with ",{aStart,bStart}];

a = aStart; fa = f[a]; dfa = df[a];
b = bStart; fb = f[b]; dfb = df[b];
fzero = f[0]; dfzero = df[0];

"We now check Fletcher's conditions";

If[fa > fzero + a * rho *dfzero,
Print["Initial bracket doesn't obey Fletcher's a value condition (i)"]];

If[(b -a )*dfa < 0 && Abs[dfa] > -sigma*dfzero,a=a,
	Print["Initial bracket doesn't obey Fletcher's derivative condition (ii)"]];

If[fb > fzero + b*rho*dfzero || fb ≥ fa, a = a,
	Print["Initial bracket doesn't obey Fletcher's b value condition (iii)"];
	Print["b=",b,"fb = ",fb," fzero + b*rho*dfzero =",fzero + b*rho*dfzero ];
	Print["fb = ",fb," fa = ",fa];
];

"These values could have been passed from FletcherMin, so we don't count them";

For[j=1,j<100,j++,

"Now we choose a new intermediate point in a smaller interval";
{targetA,targetB} = {a + tau2(b-a),b - tau3(b-a)};
 alpha = CubicMin[{{a,fa,dfa},{b,fb,dfb}}];
Prinp[alpha," ",a," ",b];
 alpha = If[alpha > targetB,targetB,alpha];
 alpha = If[alpha < targetA,targetA,alpha];
Prinp[alpha];
 falpha = f[alpha]; dfalpha = df[alpha]; 
 FEvals++; dFEvals++;

If[falpha > fzero + rho * alpha*dfzero || falpha ≥ fa,

	"alpha is still too far out. We use it as the new endpoint.";
	{b,fb,dfb} = {alpha,falpha,dfalpha};
,
	If[Abs[dfalpha] ≤ -sigma * dfzero,Return[{alpha,falpha}]];
	
	"If not, we update the interval to (alpha,a) or (alpha,b)."; 
	If[(b -a)*dfalpha ≥ 0, {b,fb,dfb} = {a,fa,dfa}];
	{a,fa,dfa} = {alpha,falpha,dfalpha};

];

];

Print["Too many iterations in FletcherBracket"];

];

GoldsteinConditions[f_,df_,a_] := Module[{rho = 0.01,sigma = 0.9},
Print["Checking modified Goldstein conditions for a = ",a];

Print ["   Abs[df[a]] < -sigma * df[0] : ", Abs[df[a]], " < ", -sigma*df[0]];
Print["   f[a] <= f[0] + a * rho * df[0] : ",f[a], " < ", f[0] + a*rho*df[0]];

TrueQ[Abs[df[a]] < -sigma*df[0] && f[a] < f[0] + a *rho*df[0]]
]

Experiments with the Fletcher Algorithm

We now test the Fletcher algorithm on some standard test cases. First, consider Sin:

Plot[Sin[x],{x,3,7}]
Sin[3.0]

f[x_] := Sin[x-3];
df[x_] := Cos[x-3];

FletcherMin[f,df,-0.95]
FEvals + dFEvals

GoldsteinConditions[f,df,1]