https://github.com/kul-forbes/ForBES
Revision dc65c2366e10a76b7166228adb65acb9a20d68a0 authored by Pantelis Sopasakis on 05 March 2018, 20:56:03 UTC, committed by GitHub on 05 March 2018, 20:56:03 UTC
crash of Forbes on Matlab 2017a
Tip revision: dc65c2366e10a76b7166228adb65acb9a20d68a0 authored by Pantelis Sopasakis on 05 March 2018, 20:56:03 UTC
Merge pull request #3 from Zilleplus/master
Merge pull request #3 from Zilleplus/master
Tip revision: dc65c23
LineSearch_MoreThuente.m
function [cachet, t, cnt, info ] = LineSearch_MoreThuente(prob, gam, cache, df0, lsopt)
% Function MoreThuenteLS
%
% The purpose of MoreThuenteLS is to find a step which satisfies
% a sufficient decrease condition and a curvature condition.
%
% At each stage MoreThuenteLS updates an interval of
% uncertainty with endpoints stx and sty. The interval of
% uncertainty is initially chosen so that it contains a
% minimizer of the modified function
%
% f(x+stp*s) - f(x) - delta*stp*(gradf(x)'s).
%
% If a step is obtained for which the modified function
% has a nonpositive function value and nonnegative derivative,
% then the interval of uncertainty is chosen so that it
% contains a minimizer of f(x+stp*s).
%
% The algorithm is designed to find a step which satisfies
% the sufficient decrease condition
%
% f(x+stp*s) <= f(x) + delta*stp*(gradf(x)'s),
%
% and the curvature condition
%
% abs(gradf(x+stp*s)'s)) <= sigma*abs(gradf(x)'s).
%
% If delta is less than sigma and if, for example, the function
% is bounded below, then there is always a step which satisfies
% both conditions. If no step can be found which satisfies both
% conditions, then the algorithm usually stops when rounding
% errors prevent further progress. In this case stp only
% satisfies the sufficient decrease condition.
%
%
% stp is a nonnegative variable. On input stp contains an
% initial estimate of a satisfactory step. On output
% stp contains the final estimate.
%
% delta and sigma are nonnegative input variables. Termination
% occurs when the sufficient decrease condition and the
% directional derivative condition are satisfied.
%
% progTol is a nonnegative input variable. Termination occurs
% when the relative width of the interval of uncertainty
% is at most progTol.
%
% tmin and tmax are nonnegative input variables which
% specify lower and upper bounds for the step.
%
% maxfev is a positive integer input variable. Termination
% occurs when the number of calls to fcn is at least
% maxfev by the end of an iteration.
%
% info is an integer output variable set as follows:
%
% info = 0 The sufficient decrease condition and the
% directional derivative condition hold.
%
% info = -2 Relative width of the interval of uncertainty
% is at most progTol.
%
% info = -1 Number of calls to fcn has reached maxfev.
%
% info = 4 The step is at the lower bound tmin.
%
% info = 5 The step is at the upper bound tmax.
%
% info = 6 Rounding errors prevent further progress.
% There may not be a step which satisfies the
% sufficient decrease and curvature conditions.
% Tolerances may be too small.
%
% nf is an integer output variable set to the number of
% calls to fcn.
%
% wa is a work array of length n.
%
% Subprograms called
%
% user-supplied......fcn
%
% MINPACK-supplied...cstep
%
% FORTRAN-supplied...abs,max,min
%
% Argonne National Laboratory. MINPACK Project. June 1983
% Jorge J. More', David J. Thuente
%
% **********
% Q, A, C, f2, proxg
cnt = [0, 0, 0, 0, 0];
nf = 0;
p5 = .5;
p66 = .66;
xtrapf = 4;
info = 0;
infoc = 1;
% Compute the initial gradient in the search direction
% and check that s is a descent direction.
%
%
% Initialize local variables.
%
t = lsopt.tau0;
f = cache.FBE;
brackt = 0;
stage1 = 1;
finit = f;
dgtest = lsopt.delta*df0;
width = lsopt.tmax - lsopt.tmin;
width1 = 2*width;
%
% The variables stx, fx, dgx contain the values of the step,
% function, and directional derivative at the best step.
% The variables sty, fy, dgy contain the value of the step,
% function, and derivative at the other endpoint of
% the interval of uncertainty.
% The variables stp, f, dg contain the values of the step,
% function, and derivative at the current step.
%
stx = 0;
fx = finit;
dgx = df0;
sty = 0;
fy = finit;
dgy = df0;
%
% Start of iteration.
%
while (1)
%
% Set the minimum and maximum steps to correspond
% to the present interval of uncertainty.
%
if (brackt)
stmin = min(stx,sty);
stmax = max(stx,sty);
else
stmin = stx;
stmax = t + xtrapf*(t - stx);
end
%
% Force the step to be within the bounds tmax and tmin.
%
t = max(t,lsopt.tmin);
t = min(t,lsopt.tmax);
%
% If an unusual termination is to occur then let
% t be the lowest point obtained so far.
%
if ((brackt && (t <= stmin || t >= stmax)) || nf >= lsopt.maxfev-1 || infoc == 0 || (brackt && stmax-stmin <= lsopt.progTol*stmax))
t = stx;
end
%
% Evaluate the function and gradient at t
% and compute the directional derivative.
%
[cachet, cnt1] = LineFBE(prob, gam, t, cache, 3);
cnt = cnt+cnt1;
nf = nf+1;
f = cachet.FBE; dg = cachet.dFBE;
ftest1 = finit + t*dgtest;
%
% Test for convergence.
%
if ((brackt && (t <= stmin || t >= stmax)) || infoc == 0)
% Rounding errors prevent further progress
info = 6;
return
end
if (t == lsopt.tmax && f <= ftest1 && dg <= dgtest)
% The step is at the upper bound tmax
info = 5;
return
end
if (t == lsopt.tmin && (f > ftest1 || dg >= dgtest))
% The step is at the lower bound tmin
info = 4;
return
end
if (nf >= lsopt.maxfev)
info = -1;
return
end
if (brackt && stmax-stmin <= lsopt.progTol*stmax)
info = -2;
return
end
if (f <= ftest1 && abs(dg) <= lsopt.sigma*(-df0))
info = 0;
return
end
%
% In the first stage we seek a step for which the modified
% function has a nonpositive value and nonnegative derivative.
%
if (stage1 && f <= ftest1 && dg >= min(lsopt.delta,lsopt.sigma)*df0)
stage1 = 0;
end
%
% A modified function is used to predict the step only if
% we have not obtained a step for which the modified
% function has a nonpositive function value and nonnegative
% derivative, and if a lower function value has been
% obtained but the decrease is not sufficient.
%
if (stage1 && f <= fx && f > ftest1)
%
% Define the modified function and derivative values.
%
fm = f - t*dgtest;
fxm = fx - stx*dgtest;
fym = fy - sty*dgtest;
dgm = dg - dgtest;
dgxm = dgx - dgtest;
dgym = dgy - dgtest;
%
% Call cstep to update the interval of uncertainty
% and to compute the new step.
%
[stx,fxm,dgxm,sty,fym,dgym,t,fm,dgm,brackt,infoc] = MoreThuenteCstep(stx,fxm,dgxm,sty,fym,dgym,t,fm,dgm,brackt,stmin,stmax);
%
% Reset the function and gradient values for f.
%
fx = fxm + stx*dgtest;
fy = fym + sty*dgtest;
dgx = dgxm + dgtest;
dgy = dgym + dgtest;
else
%
% Call cstep to update the interval of uncertainty
% and to compute the new step.
%
[stx,fx,dgx,sty,fy,dgy,t,f,dg,brackt,infoc] = MoreThuenteCstep(stx,fx,dgx,sty,fy,dgy,t,f,dg,brackt,stmin,stmax);
end
%
% Force a sufficient decrease in the size of the
% interval of uncertainty.
%
if (brackt)
if (abs(sty-stx) >= p66*width1)
t = stx + p5*(sty - stx);
end
width1 = width;
width = abs(sty-stx);
end
%
% End of iteration.
%
end
end
function [stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,info] = MoreThuenteCstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,stpmin,stpmax)
% Subroutine cstep
%
% The purpose of cstep is to compute a safeguarded step for
% a linesearch and to update an interval of uncertainty for
% a minimizer of the function.
%
% The parameter stx contains the step with the least function
% value. The parameter stp contains the current step. It is
% assumed that the derivative at stx is negative in the
% direction of the step. If brackt is set true then a
% minimizer has been bracketed in an interval of uncertainty
% with endpoints stx and sty.
%
% The subroutine statement is
%
% subroutine cstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt,
% stpmin,stpmax,info)
%
% where
%
% stx, fx, and dx are variables which specify the step,
% the function, and the derivative at the best step obtained
% so far. The derivative must be negative in the direction
% of the step, that is, dx and stp-stx must have opposite
% signs. On output these parameters are updated appropriately.
%
% sty, fy, and dy are variables which specify the step,
% the function, and the derivative at the other endpoint of
% the interval of uncertainty. On output these parameters are
% updated appropriately.
%
% stp, fp, and dp are variables which specify the step,
% the function, and the derivative at the current step.
% If brackt is set true then on input stp must be
% between stx and sty. On output stp is set to the new step.
%
% brackt is a logical variable which specifies if a minimizer
% has been bracketed. If the minimizer has not been bracketed
% then on input brackt must be set false. If the minimizer
% is bracketed then on output brackt is set true.
%
% stpmin and stpmax are input variables which specify lower
% and upper bounds for the step.
%
% info is an integer output variable set as follows:
% If info = 1,2,3,4,5, then the step has been computed
% according to one of the five cases below. Otherwise
% info = 0, and this indicates improper input parameters.
%
% Subprograms called
%
% FORTRAN-supplied ... abs,max,min,sqrt
% ... dble
%
% Argonne National Laboratory. MINPACK Project. June 1983
% Jorge J. More', David J. Thuente
%
% **********
p66 = 0.66;
info = 0;
%
% Check the input parameters for errors.
%
if ((brackt && (stp <= min(stx,sty) || stp >= max(stx,sty))) || dx*(stp-stx) >= 0.0 || stpmax < stpmin)
return
end
%
% Determine if the derivatives have opposite sign.
%
sgnd = dp*(dx/abs(dx));
%
% First case. A higher function value.
% The minimum is bracketed. If the cubic step is closer
% to stx than the quadratic step, the cubic step is taken,
% else the average of the cubic and quadratic steps is taken.
%
if (fp > fx)
info = 1;
bound = 1;
theta = 3*(fx - fp)/(stp - stx) + dx + dp;
s = norm([theta,dx,dp],inf);
gamma = s*sqrt((theta/s)^2 - (dx/s)*(dp/s));
if (stp < stx)
gamma = -gamma;
end
p = (gamma - dx) + theta;
q = ((gamma - dx) + gamma) + dp;
r = p/q;
stpc = stx + r*(stp - stx);
stpq = stx + ((dx/((fx-fp)/(stp-stx)+dx))/2)*(stp - stx);
if (abs(stpc-stx) < abs(stpq-stx))
stpf = stpc;
else
stpf = stpc + (stpq - stpc)/2;
end
brackt = 1;
%
% Second case. A lower function value and derivatives of
% opposite sign. The minimum is bracketed. If the cubic
% step is closer to stx than the quadratic (secant) step,
% the cubic step is taken, else the quadratic step is taken.
%
elseif (sgnd < 0.0)
info = 2;
bound = 0;
theta = 3*(fx - fp)/(stp - stx) + dx + dp;
s = norm([theta,dx,dp],inf);
gamma = s*sqrt((theta/s)^2 - (dx/s)*(dp/s));
if (stp > stx)
gamma = -gamma;
end
p = (gamma - dp) + theta;
q = ((gamma - dp) + gamma) + dx;
r = p/q;
stpc = stp + r*(stx - stp);
stpq = stp + (dp/(dp-dx))*(stx - stp);
if (abs(stpc-stp) > abs(stpq-stp))
stpf = stpc;
else
stpf = stpq;
end
brackt = 1;
%
% Third case. A lower function value, derivatives of the
% same sign, and the magnitude of the derivative decreases.
% The cubic step is only used if the cubic tends to infinity
% in the direction of the step or if the minimum of the cubic
% is beyond stp. Otherwise the cubic step is defined to be
% either stpmin or stpmax. The quadratic (secant) step is also
% computed and if the minimum is bracketed then the the step
% closest to stx is taken, else the step farthest away is taken.
%
elseif (abs(dp) < abs(dx))
info = 3;
bound = 1;
theta = 3*(fx - fp)/(stp - stx) + dx + dp;
s = norm([theta,dx,dp],inf);
%
% The case gamma = 0 only arises if the cubic does not tend
% to infinity in the direction of the step.
%
gamma = s*sqrt(max(0.,(theta/s)^2 - (dx/s)*(dp/s)));
if (stp > stx)
gamma = -gamma;
end
p = (gamma - dp) + theta;
q = (gamma + (dx - dp)) + gamma;
r = p/q;
if (r < 0.0 && gamma ~= 0.0)
stpc = stp + r*(stx - stp);
elseif (stp > stx)
stpc = stpmax;
else
stpc = stpmin;
end
stpq = stp + (dp/(dp-dx))*(stx - stp);
if (brackt)
if (abs(stp-stpc) < abs(stp-stpq))
stpf = stpc;
else
stpf = stpq;
end
else
if (abs(stp-stpc) > abs(stp-stpq))
stpf = stpc;
else
stpf = stpq;
end
end
%
% Fourth case. A lower function value, derivatives of the
% same sign, and the magnitude of the derivative does
% not decrease. If the minimum is not bracketed, the step
% is either stpmin or stpmax, else the cubic step is taken.
%
else
info = 4;
bound = 0;
if (brackt)
theta = 3*(fp - fy)/(sty - stp) + dy + dp;
s = norm([theta,dy,dp],inf);
gamma = s*sqrt((theta/s)^2 - (dy/s)*(dp/s));
if (stp > sty)
gamma = -gamma;
end
p = (gamma - dp) + theta;
q = ((gamma - dp) + gamma) + dy;
r = p/q;
stpc = stp + r*(sty - stp);
stpf = stpc;
elseif (stp > stx)
stpf = stpmax;
else
stpf = stpmin;
end
end
%
% Update the interval of uncertainty. This update does not
% depend on the new step or the case analysis above.
%
if (fp > fx)
sty = stp;
fy = fp;
dy = dp;
else
if (sgnd < 0.0)
sty = stx;
fy = fx;
dy = dx;
end
stx = stp;
fx = fp;
dx = dp;
end
%
% Compute the new step and safeguard it.
%
stpf = min(stpmax,stpf);
stpf = max(stpmin,stpf);
stp = stpf;
if (brackt && bound)
if (sty > stx)
stp = min(stx+p66*(sty-stx),stp);
else
stp = max(stx+p66*(sty-stx),stp);
end
end
end
Computing file changes ...