Home > source > @penlab > ls_pennon.m

# ls_pennon

## PURPOSE

rewritten line_search() @ line_search_els.c @ Pennlp v.2.3 & Pennon v.0.9

## SYNOPSIS

function [rRelStep, nFlag]=ls_pennon(obj, dir)

## DESCRIPTION

     rewritten line_search() @ line_search_els.c @ Pennlp v.2.3 & Pennon v.0.9

     changes in 'obj': xall, ALx, ALdx

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SOURCE CODE

    0001 function [rRelStep, nFlag]=ls_pennon(obj, dir)
    0002 % rewritten line_search() @ line_search_els.c @ Pennlp v.2.3 & Pennon v.0.9
    0003 %
    0004 % changes in 'obj': xall, ALx, ALdx
    0005
    0006 %%%%%%%%% Settings %%%%%%%%%%
    0007   XTOL = 1.0E-1;
    0008   ETA = 0.9;
    0009 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0010
    0011 % d = dir ... computed search direction
    0012
    0013   MAX_LS_ITER = obj.allopts.max_ls_iter;
    0014
    0015   xall0=obj.xall;   % original point
    0016   % whenever I want to evaluate at a new point, xall<-x_it
    0017   fx=obj.ALx;
    0018
    0019   % double pi0, pi, mu, divs, mu1s, mu2s, eff, kappa, fact, q;
    0020   mu1 = 1.; mu2 = 1.; sml=0.5;   %%% fbest, f, f0, f1, f2,
    0021   alp0 = 0; alp1 = 0; alp2 = 0; alps = 0; alpm = 0; alpd = 0; omega1 = 0; omega2 = 0;
    0022   stpmin = 0; stpmax = 1.0E20; alp = 1.; denom = 0.;
    0023   FTOL = 1-ETA;
    0024
    0025   % int
    0026   lsiter=0;
    0027   nEff = 0; iterations = 0;
    0028   k = 0; ksh = 0; ier = 0;
    0029   bLong = 0; bSafeguard = 0;
    0030
    0031   bLinesearch = 1;
    0032
    0033   rRelStep=0.; nFlag=99;
    0034
    0035   %if(Uequal)
    0036   %  if(Uequal->dim) {
    0037   %    stpmax=1.0;
    0038   %  }
    0039
    0040
    0041   pi=dir'*obj.ALdx;
    0042   % computation of the rest of the gradient of our merit function
    0043   %if(Uequal) {
    0044     % v Src 0.9 je pouze:  (ale to je temer jiste blbe
    0045   %  pi -= rNu * rEqNorm;
    0046     % misto nasledujiciho z Pennlp 2.3(fi)
    0047   %  spm3v_mlt( (SPMAT3 *) gradeq, Cequal, v_tmp_short );
    0048   %  rTmp = in_prod(v_tmp_short, d);
    0049   %  pi += rNu*rTmp;
    0050   %  printf("pi: %e, rTmp: %e\n", pi, rTmp);
    0051   %}
    0052
    0053   %if(Uequal) {
    0054   %  v_copy(Uequal, u_it);
    0055   %  for(i=0; i<Uequal->dim; i++)
    0056   %    u_d->ve[i] = d->ve[X->dim+i];
    0057   %    // smer pro multiplikatory je "schovany" v d za x, ale dim je jen pro x
    0058   %}
    0059
    0060   %if(Uequal) {
    0061   %  spm3v_mlt( (SPMAT3 *) gradeq, u_d, v_tmp_short );
    0062   %  rTmp = in_prod(v_tmp_short, d);
    0063   %  pi += rTmp;
    0064   %  printf("pi2: %e, dAv: %e", pi, rTmp);
    0065   %}
    0066
    0067   % Is it really descent
    0068   if (pi > 0)
    0069     obj.print(3,Inf,'LS (pen): No Descent');
    0070
    0071     bLinesearch = 0;
    0072     rRelStep=0;
    0073     nFlag=1;
    0074     return;
    0075   elseif (abs(pi) < 1.0E-10)
    0076     alp = 1.0;
    0077     %obj.xall = xall0 + alp*dir;   % not needed here
    0078     %if (Uequal)
    0079     %  UEQ = u_it + alp*u_dir;
    0080     %end
    0081
    0082     %while((*fx =  (* objective) (x_it)) >= HUGE_VAL) {
    0083     %  alp = alp/2;
    0084     %  x_it=x+alp*dir;
    0085     %}
    0086
    0087     if (abs(pi) < 1.0E-15)
    0088       obj.ls_small_cnt=obj.ls_small_cnt+1;
    0089     else
    0090       obj.ls_small_cnt=0;
    0091     end
    0092
    0093     if (obj.ls_small_cnt > 2)
    0094       % x remains unchanged --> nothing todo
    0095       % but objective() could be called elsewhere then at x
    0096          % doesn't seem so at this version (JF 4/11/11)
    0097       %fx = fnc.obj(x);
    0098       rRelStep=0.;
    0099       nFlag=2;
    0100       obj.print(3,Inf,'LS (pen):  %.4e (almost not descent), no step', pi);
    0101     else        % still acceptable
    0102       obj.xall = xall0 + alp*dir;
    0103       %x=x_it;
    0104       obj.eval_alx();
    0105       obj.eval_aldx();
    0106       %fx = fnc.obj(x);
    0107       %gradx = fnc.obj_grad(x);
    0108       obj.print(3,Inf,'LS (pen):  %.4e (small full step)', pi);
    0109       rRelStep=1.;
    0110       nFlag=0;
    0111       lsiter=lsiter+1;
    0112     end
    0113
    0114     return;
    0115   %else
    0116   %  obj.print(3,Inf,'LS (pen):  %.4e', pi);   %%% ?
    0117   end
    0118
    0119   obj.ls_small_cnt=0;
    0120
    0121   %*****************************************************
    0122   %* Initialize line search ... *
    0123   %*****************************************************
    0124
    0125   if (FTOL <= 0)
    0126     mu1s = 1./3.;
    0127     mu2s = 2./3.;
    0128     eff = 4./9.;
    0129   else
    0130     mu1s = min(FTOL,1/3);
    0131     mu2s = 1-mu1s;
    0132     eff = 2*mu1s*mu2s;
    0133   end
    0134
    0135   if (XTOL >= 1.)
    0136     kappa = 2.;
    0137   else
    0138     kappa = 1. / min(1.,1.-XTOL);
    0139   end
    0140
    0141   f0 = obj.ALx; f1 = obj.ALx; f2 = obj.ALx; fbest = obj.ALx;
    0142   pi0 = pi;
    0143
    0144   %*****************************************************
    0145   %* actual line search iteration ... *
    0146   %*****************************************************
    0147
    0148   nFlag=3;   % max iter
    0149   % for equalities should be:  iterations < LS_MAX_ITER && (< LS_MAX_EQ_ITER if equal)
    0150   for iterations=1:MAX_LS_ITER
    0151     ier = 0;
    0152     % get new function value
    0153     %//printf(" * %f", alp);
    0154
    0155     %x_it = x+alp*dir;
    0156     obj.xall = xall0+alp*dir;
    0157     %if(Uequal)
    0158     %  UEQ = u_it + alp*u_dir;
    0159     %end
    0160
    0161     %f = fnc.obj(x_it);
    0162     obj.eval_alx();
    0163     f = obj.ALx;
    0164     nEff = 0;
    0165
    0166     %*****************************************************
    0167     %* get new trial step alp
    0168     %*****************************************************
    0169
    0170     % count requested function evaluations
    0171     k=k+1;
    0172
    0173     % check correctness of step
    0174     % note that alp1>=alp2 unless alp1=0
    0175     if (alp < stpmin || alp <= alp2)
    0176       % new step size too small
    0177       obj.print(4,Inf,'LS (pen): ERROR: new step size too small');
    0178       nEff = 0;
    0179       ier = 1;
    0180       break;
    0181     elseif (alp > stpmax || (alp >= alp1 && alp1 > 0))
    0182       % new step size too large
    0183       obj.print(4,Inf,'LS (pen): ERROR: new step size too large');
    0184       nEff = 0;
    0185       ier = 1;
    0186       break;
    0187     end
    0188
    0189     % update best point
    0190     if (f < fbest || (f == fbest && f == f0))
    0191       fbest=f;
    0192       if (f == f0 && (ksh == 0 || alp < alps))
    0193         % store smallest step with f=f0 for possible restauration
    0194         alps = alp;
    0195       end
    0196       if (f == f0 && (ksh == 0 || alp > alpm))
    0197         % store smallest step with f=f0 for possible restauration
    0198         alpm = alp;
    0199       end
    0200       if (f == f0)
    0201         ksh=ksh+1;
    0202       end
    0203     end
    0204
    0205     % get parabolic minimizer alp0
    0206     if (alp1*alp2 > 0)
    0207       alp0 = alp1;    % default if minimizer does not exist
    0208       denom = (f-f1)/(alp-alp1)-(f-f2)/(alp-alp2);
    0209       if (denom > 0)
    0210         alp0 = 0.5*(alp1+alp2-(f1-f2)/denom);
    0211       end
    0212     end
    0213
    0214     % compute Goldstein quotient
    0215     mu = (f-f0)/(alp*pi0);
    0216     if (f == f0 && f2 == f0)
    0217       % check whether f=f0 due to very short step
    0218       if (f1 < f0 || ksh <= 6)
    0219         % assume very short step
    0220         mu = 1.;
    0221       end
    0222     end
    0223
    0224     % update bracket, with safeguard when f=f0
    0225     bLong = (mu < mu1s); %%% ??? || f >= HUGE_VAL/*_isnan(f)*/);
    0226     if (bLong)
    0227       % long step, overwrite alp1
    0228       alp1 = alp;
    0229       f1 = f;
    0230       mu1 = mu;
    0231       if (k == 1 && mu > 0.5)
    0232         mu1 = 0;
    0233       end
    0234       omega1 = 2*abs(1.-mu1)/alp1;
    0235     else
    0236       % good step or short step, overwrite alp2
    0237       alp2 = alp;
    0238       f2 = f;
    0239       mu2 = mu;
    0240       omega2 = 2*abs(1.-mu2)/alp2;
    0241     end
    0242
    0243     if (f2 == f0 && f1 >= f0 && k >= 2)
    0244       % check whether f=f0 steps should be considered long
    0245       % ksh=7 indicates f=f0 for all large alp
    0246       % f1=f0 && alp1=alp2*2 indicates f2=f0 by chance
    0247       if (ksh == 7 || (ksh == 1 && f1 > f0 && alp1 == alp2*1.1))
    0248         ksh=9999;      % ends f=f0 tests
    0249         % restore alp1=alps, alp2=0
    0250         alp1 = alps;
    0251         f1 = f0;
    0252         mu1 = 0;
    0253         omega1 = 0;
    0254         alp2 = 0;
    0255         f2 = f0;
    0256         mu2 = 1;
    0257         omega2 = 0;
    0258       end
    0259     end
    0260
    0261     %**********************************
    0262     %* convergence tests              *
    0263     %**********************************
    0264
    0265     %if (f0 >= HUGE_VAL && f < f0)
    0266       % f0=infinity and f is finite
    0267     %  obj.print(4,Inf,'WARNING: infinite f0 improved to finite');
    0268     %  nEff = 2;
    0269     %  ier = 0;
    0270     %  break;
    0271     %end
    0272
    0273     if (bLong)
    0274       if (alp1 == stpmin && alp1 ~= 0.)
    0275         % stpmin too large
    0276         obj.print(3,Inf,'WARNING: end of line search forced by stpmin');
    0277         nEff = 2;
    0278         ier = 0;
    0279         break;
    0280       end
    0281     else
    0282       % now the sufficient decrease condition holds
    0283       if (0)   %/*_isnan(f1)*/
    0284         if (k >= 3 && alp == alp2 && alp1 <= kappa*alp2)
    0285           % step close to a NaN point !!!!!!
    0286           nEff = 1;
    0287           ier = 0;
    0288           break;
    0289         end
    0290       end
    0291
    0292       if (alp2*mu2*max(omega1,omega2) >= eff && f ~= f0)
    0293         % efficiency condition holds
    0294         nEff = 1;
    0295         ier = 0;
    0296         break;
    0297       end
    0298
    0299       % check step size restrictions
    0300       if (alp1 == stpmin && alp1 ~= 0.)
    0301         % stpmin too large
    0302         obj.print(3,Inf,'WARNING: end of line search forced by stpmin');
    0303         nEff = 2;
    0304         ier = 0;
    0305         break;
    0306       end
    0307
    0308       if (alp2 == stpmax)
    0309         % stpmax too small
    0310         obj.print(3,Inf,'WARNING: end of line search forced by stpmax');
    0311         nEff = 2;
    0312         ier = 0;
    0313         break;
    0314       end
    0315     end
    0316
    0317     %********************************************
    0318     %* compute new trial stepsize (goto150)     *
    0319     %********************************************
    0320
    0321     bSafeguard = 1;
    0322
    0323     % extrapolation phase
    0324     if (alp1 == 0.)
    0325       % only short steps, mu>1/2, expand
    0326       if (fbest == f0)
    0327         % escape roundoff
    0328         alp = alp2/sml;
    0329         % safeguard in case fbest=f0 by chance
    0330         if (ksh == 1)
    0331           alp = alp2*1.1;
    0332         end
    0333       elseif (k == 1)
    0334         % second trial point, safeguarded parabolic minimization
    0335         alp = alp/max(2-2*mu,sml);
    0336       else
    0337         % expand by divs=cut(4*mu2,[2,4])
    0338         divs = max(2.,min(4.,4.*mu2));
    0339         alp = alp2*divs;
    0340       end
    0341
    0342       bSafeguard = 0;
    0343     else
    0344       % interpolation phase
    0345       if (alp2 == 0.)
    0346         % only long steps, mu<1/2, contract
    0347         % linear interpolation to mu=0.5
    0348         alpd = alp1/(2-2*mu1);
    0349       elseif (fbest == f0)
    0350         % escape roundoff by geometric mean step
    0351         alp = alp2*sqrt(alp1/alp2);
    0352         % safeguard in case fbest=f0 by chance
    0353         if (ksh == 1)
    0354           alp = min(alp,alp2*1.1);
    0355         end
    0356         bSafeguard = 0;
    0357       else
    0358         % test interiority condition for parabolic minimizer alp0
    0359         fact = (alp1/alp2)^(1/3);
    0360         if (denom > 0 && (alp0 >= fact*alp2 && alp1 >= fact*alp0))
    0361           % safeguarded parabolic minimization
    0362           alpd = alp0;
    0363         else
    0364           alpd = alp2;    % is updated by safeguards
    0365         end
    0366       end
    0367     end
    0368
    0369     if (bSafeguard)
    0370       % compute relative step location
    0371       q = (alpd-alp2)/(alp1-alp2);
    0372
    0373       %*********************************************
    0374       %* safeguards                                *
    0375       %*********************************************
    0376
    0377       if (alp1 == 0)
    0378         error('programming error in extrapolation step of line search (alp1==0)');
    0379       elseif (alp2 == 0)
    0380         % only long steps mu<1/2
    0381         if (q >= sml*sml && q <= sml)
    0382           % contract by factor q <= sml
    0383           alp = q*alp1;
    0384         else
    0385           % contract by sml
    0386           alp = sml*alp1;
    0387         end
    0388       else
    0389         % force interiority condition for interpolation step
    0390         alpd = alp2+q*(alp1-alp2);
    0391         fact = (alp1/alp2)^(1/3);
    0392         if (0)      %/*_isnan(f1)*/
    0393           % f1=NaN, take a short interior step
    0394           alpd = fact*alp2;
    0395         else
    0396           % f1 is a valid number or infinite
    0397           if (alpd<fact*alp2)
    0398             % cut interpolation step
    0399             alpd = fact*alp2;
    0400           end
    0401         end
    0402
    0403         alp = min(alpd,alp1/fact);
    0404         if (alp ~= alpd || alp >= alp1 || alp<=alp2)
    0405           % step not interior enough, use geometric mean step
    0406           alp = alp1*sqrt(alp2/alp1);
    0407         end
    0408       end
    0409     end     % of if(bSafeguard)
    0410
    0411     %*********************************************
    0412     %* check bounds and assign task  (goto2000)  *
    0413     %*********************************************
    0414
    0415     % check bounds on step size
    0416     if (alp < stpmin)
    0417       alp = stpmin;
    0418     elseif (alp > stpmax)
    0419       alp=stpmax;
    0420     end
    0421
    0422     if (alp > alp2 && (alp < alp1 || alp1 == 0.))
    0423       % good trial step; request next function value
    0424       % task='VALUE'
    0425       ;
    0426     else
    0427       % interval too small for splitting
    0428       if (alp2*mu2*max(omega1,omega2) >= eff && f ~= f0)
    0429         % efficiency condition holds
    0430         obj.print(3,Inf,'WARNING: end of line search forced by machine accuracy');
    0431         nEff = 2;
    0432         ier = 0;
    0433         break;
    0434       else
    0435         obj.print(3,Inf,'WARNING: end of line search forced by machine accuracy');
    0436         nEff = 0;
    0437         ier = 1;
    0438         break;
    0439       end
    0440     end
    0441
    0442     if (alp2 == 0)
    0443       fact = alp1/alp;
    0444     else
    0445       fact = alp/alp2;
    0446     end
    0447
    0448     %******************************************************
    0449
    0450   end    % for iterations
    0451
    0452
    0453   %*****************************************************
    0454
    0455   if (ier == 0)
    0456     if (f < fx)   %%% by mela byt merit??
    0457       % x=x_it;   % already done JF041111
    0458       bLinesearch = 0;
    0459       % fx = fnc.obj(x);
    0460       fx=f;
    0461       obj.eval_aldx();
    0462       %gradx = fnc.obj_grad(x);
    0463       nFlag=0;
    0464     else
    0465       if (ksh>0 && alpm > 0.) %fbest < HUGE_VAL && alpm > 0.)
    0466         % We have descent direction, but
    0467         % function values cannot be improved due to
    0468         % numerical inexactness
    0469         alp = alpm;
    0470         %x_it=x+alp*dir;
    0471         %x=x_it;
    0472         obj.xall=xall0+alp*dir;
    0473         %if (Uequal)
    0474     %  UEQ = u_it + alp*u_dir;
    0475         %end
    0476
    0477         bLinesearch = 0;
    0478
    0479         %fx = fnc.obj(x);
    0480         %gradx = fnc.obj_grad(x);
    0481         obj.eval_alx();
    0482         obj.eval_aldx();
    0483         fx = obj.ALx;
    0484         %gradx = fnc.obj_grad(x);
    0485     nFlag=0;
    0486       else
    0487         alp = 0;
    0488       end
    0489     end
    0490   else
    0491     alp = 0;
    0492   end
    0493
    0494   lsiter=lsiter+iterations;
    0495   obj.lsiter_last=obj.lsiter_last + lsiter;
    0496
    0497   obj.print(3,Inf,'LS (pen): %.4e, %d steps, Step width: %f', pi0, iterations, alp);
    0498
    0499   %return alp;
    0500   rRelStep=alp;
    0501   return;
    0502
    0503