stogo/local.cc

   1
   2 /*
   3    Local search - A trust region algorithm with BFGS update.
   4 */
   5
   6 #include <iostream>
   7 #include <stdlib.h>
   8
   9 #include "stogo_config.h"
  10 #include "global.h"
  11 #include "local.h"
  12 #include "tools.h"
  13
  14 #ifdef NLOPT_UTIL_H
  15 #  define IF_NLOPT_CHECK_EVALS stop->nevals++; \
  16                                if (nlopt_stop_evalstime(stop)) \
  17                                   return LS_MaxEvalTime
  18 #else
  19 #  define IF_NLOPT_CHECK_EVALS
  20 #endif
  21
  22 ////////////////////////////////////////////////////////////////////////
  23 // SGJ, 2007: allow local to use local optimizers in NLopt, to compare
  24 // to the BFGS code below
  25 #if 0
  26 #include "nlopt.h"
  27
  28 typedef struct {
  29   Global *glob;
  30   double maxgrad;
  31   nlopt_stopping *stop;
  32 } f_local_data;
  33
  34 static double f_local(int n, const double *x, double *grad, void *data_)
  35 {
  36   f_local_data *data = (f_local_data *) data_;
  37   double f;
  38   RVector xv, gv;
  39   // hack to avoid pointless copy of x and grad
  40   xv.len = gv.len = n;
  41   xv.elements = const_cast<double *>(x);
  42   gv.elements = grad;
  43   f=data->glob->ObjectiveGradient(xv, gv,
  44                                    grad?OBJECTIVE_AND_GRADIENT:OBJECTIVE_ONLY);
  45   if (grad) data->maxgrad = max(data->maxgrad, normInf(gv));
  46   xv.elements = gv.elements = 0; // prevent deallocation
  47   data->stop->nevals++;
  48   return f;
  49 }
  50 #endif
  51 ////////////////////////////////////////////////////////////////////////
  52
  53 int local(Trial &T, TBox &box, TBox &domain, double eps_cl, double *mgr,
  54           Global &glob, int axis, RCRVector x_av
  55 #ifdef NLOPT_UTIL_H
  56           , nlopt_stopping *stop
  57 #endif
  58           ) {
  59
  60   int n=box.GetDim();
  61   RVector x(n);
  62   double tmp, f;
  63
  64   x=T.xvals ;
  65
  66 #ifdef LS_DEBUG
  67   cout << "Local Search, x=" << x << endl;
  68 #endif
  69
  70   if (box.OutsideBox(x, domain) != 0) {
  71     cout << "Starting point is not inside the boundary. Exiting...\n" ;
  72     exit(1) ;
  73     return LS_Out ;
  74   }
  75
  76   // Check if we are close to a stationary point located previously
  77   if (box.CloseToMin(x, &tmp, eps_cl)) {
  78 #ifdef LS_DEBUG
  79      cout << "Close to a previously located stationary point, exiting" << endl;
  80 #endif
  81      T.objval=tmp;
  82      return LS_Old ;
  83    }
  84
  85 #if 0
  86
  87   if (axis != -1) {
  88     cout << "NLopt code only works with axis == -1, exiting...\n" ;
  89     exit(EXIT_FAILURE);
  90   }
  91   f_local_data data;
  92   data.glob = &glob;
  93   data.maxgrad = *mgr;
  94   data.stop = stop;
  95   nlopt_result ret = nlopt_minimize(NLOPT_LOCAL_LBFGS, n, f_local, &data,
  96                                     box.lb.raw_data(), box.ub.raw_data(),
  97                                     x.raw_data(), &f,
  98                                     stop->minf_max,
  99                                     stop->ftol_rel, stop->ftol_abs,
 100                                     stop->xtol_rel, stop->xtol_abs,
 101                                     stop->maxeval - stop->nevals,
 102                                     stop->maxtime - stop->start);
 103   *mgr = data.maxgrad;
 104   T.xvals=x ; T.objval=f ;
 105   if (ret == NLOPT_MAXEVAL_REACHED || ret == NLOPT_MAXTIME_REACHED)
 106     return LS_MaxEvalTime;
 107   else if (ret > 0)
 108     return LS_New;
 109   else
 110     return LS_Out; // failure
 111
 112 #else /* not using NLopt local optimizer ... use original STOgo BFGS code */
 113
 114   int k_max, info, outside = 0;
 115   int k, i, good_enough, iTmp ;
 116
 117   double maxgrad, delta, f_new;
 118   double alpha, gamma, beta, d2, s2, nom, den, ro ;
 119   double nrm_sd, nrm_hn, snrm_hn, nrm_dl ;
 120   RVector g(n), h_sd(n), h_dl(n), h_n(n), x_new(n), g_new(n) ;
 121   RVector s(n),y(n),z(n),w(n) ; // Temporary vectors
 122   RMatrix B(n), H(n) ;          // Hessian and it's inverse
 123
 124   k_max = max_iter*n ;
 125
 126   // Initially B and H are equal to the identity matrix
 127   B=0 ; H=0 ;
 128   for (i=0 ; i<n ; i++) {
 129     B(i,i)=1 ;
 130     H(i,i)=1 ;
 131   }
 132
 133   RVector g_av(x_av.GetLength());
 134   if (axis==-1) {
 135     f=glob.ObjectiveGradient(x,g,OBJECTIVE_AND_GRADIENT);
 136   }
 137   else {
 138     x_av(axis)=x(0);
 139     f=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
 140     g(0)=g_av(axis);
 141   }
 142   IF_NLOPT_CHECK_EVALS;
 143   FC++;GC++;
 144
 145   if (axis == -1) {
 146     // Skipping AV
 147 #ifdef INI3
 148     // Elaborate scheme to initalize delta
 149     delta=delta_coef*norm2(g) ;
 150     copy(g,z) ;
 151     axpy(1.0,x,z) ;
 152     if (!box.InsideBox(z)) {
 153       if (box.Intersection(x,g,z)==TRUE) {
 154         axpy(-1.0,x,z) ;
 155         delta=min(delta,delta_coef*norm2(z)) ;
 156       }
 157       else {
 158         // Algorithm broke down, use INI1
 159         delta = (1.0/7)*box.ShortestSide(&iTmp) ;
 160       }
 161     }
 162 #endif
 163 #ifdef INI2
 164     // Use INI2 scheme
 165     delta = box.ClosestSide(x)*delta_coef ;
 166     if (delta<MacEpsilon)
 167       // Patch to avoid trust region with radius close to zero
 168       delta = (1.0/7)*box.ShortestSide(&iTmp) ;
 169 #endif
 170 #ifdef INI1
 171     delta = delta_coef*box.ShortestSide(&iTmp) ;
 172 #endif
 173   }
 174   else {
 175     // Use a simple scheme for the 1D minimization (INI1)
 176     delta = (1.0/7.0)*box.ShortestSide(&iTmp) ;
 177   }
 178
 179   k=0 ; good_enough = 0 ; info=LS_New ; outside=0 ;
 180   maxgrad=*mgr ;
 181   while (good_enough == 0) {
 182     k++ ;
 183     if (k>k_max) {
 184 #ifdef LS_DEBUG
 185       cout << "Maximum number of iterations reached\n" ;
 186 #endif
 187       info=LS_MaxIter ;
 188       break ;
 189     }
 190
 191     // Update maximal gradient value
 192     maxgrad=max(maxgrad,normInf(g)) ;
 193
 194     // Steepest descent, h_sd = -g
 195     copy(g,h_sd) ;
 196     scal(-1.0,h_sd) ;
 197     nrm_sd=norm2(h_sd) ;
 198
 199     if (nrm_sd < epsilon) {
 200       // Stop criterion (gradient) fullfilled
 201 #ifdef LS_DEBUG
 202       cout << "Gradient small enough" << endl ;
 203 #endif
 204       good_enough = 1 ;
 205       break ;
 206     }
 207
 208     // Compute Newton step, h_n = -H*g
 209     gemv('N',-1.0, H, g, 0.0, h_n) ;
 210     nrm_hn = norm2(h_n) ;
 211
 212     if (nrm_hn < delta) {
 213       // Pure Newton step
 214       copy(h_n, h_dl) ;
 215 #ifdef LS_DEBUG
 216       cout << "[Newton step]      " ;
 217 #endif
 218     }
 219     else {
 220       gemv('N',1.0,B,g,0.0,z) ;
 221       tmp=dot(g,z) ;
 222       if (tmp==0) {
 223         info = LS_Unstable ;
 224         break ;
 225       }
 226       alpha=(nrm_sd*nrm_sd)/tmp ; // Normalization (N38,eq. 3.30)
 227       scal(alpha,h_sd) ;
 228       nrm_sd=fabs(alpha)*nrm_sd ;
 229
 230       if (nrm_sd >= delta) {
 231         gamma = delta/nrm_sd ; // Normalization (N38, eq. 3.33)
 232         copy(h_sd,h_dl) ;
 233         scal(gamma,h_dl) ;
 234 #ifdef LS_DEBUG
 235         cout << "[Steepest descent]  " ;
 236 #endif
 237       }
 238       else {
 239         // Combination of Newton and SD steps
 240         d2 = delta*delta ;
 241         copy(h_sd,s) ;
 242         s2=nrm_sd*nrm_sd ;
 243         nom = d2 - s2 ;
 244         snrm_hn=nrm_hn*nrm_hn ;
 245         tmp = dot(h_n,s) ;
 246         den = tmp-s2 + sqrt((tmp-d2)*(tmp-d2)+(snrm_hn-d2)*(d2-s2)) ;
 247         if (den==0) {
 248           info = LS_Unstable ;
 249           break ;
 250         }
 251         // Normalization (N38, eq. 3.31)
 252         beta = nom/den ;
 253         copy(h_n,h_dl) ;
 254         scal(beta,h_dl) ;
 255         axpy((1-beta),h_sd,h_dl) ;
 256 #ifdef LS_DEBUG
 257         cout << "[Mixed step]        " ;
 258 #endif
 259       }
 260     }
 261     nrm_dl=norm2(h_dl) ;
 262
 263     //x_new = x+h_dl ;
 264     copy(x,x_new) ;
 265     axpy(1.0,h_dl,x_new) ;
 266
 267     // Check if x_new is inside the box
 268     iTmp=box.OutsideBox(x_new, domain) ;
 269     if (iTmp == 1) {
 270 #ifdef LS_DEBUG
 271       cout << "x_new is outside the box " << endl ;
 272 #endif
 273       outside++ ;
 274       if (outside>max_outside_steps) {
 275         // Previous point was also outside, exit
 276         break ;
 277       }
 278     }
 279     else if (iTmp == 2) {
 280 #ifdef LS_DEBUG
 281       cout << " x_new is outside the domain" << endl ;
 282 #endif
 283       info=LS_Out ;
 284       break ;
 285     }
 286     else {
 287       outside=0 ;
 288     }
 289
 290     // Compute the gain
 291     if (axis==-1)
 292       f_new=glob.ObjectiveGradient(x_new,g_new,OBJECTIVE_AND_GRADIENT);
 293     else {
 294       x_av(axis)=x_new(0);
 295       f_new=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
 296     }
 297     IF_NLOPT_CHECK_EVALS;
 298     FC++; GC++;
 299     gemv('N',0.5,B,h_dl,0.0,z);
 300     ro = (f_new-f) / (dot(g,h_dl) + dot(h_dl,z)); // Quadratic model
 301     if (ro > 0.75) {
 302       delta = delta*2;
 303     }
 304     if (ro < 0.25) {
 305       delta = delta/3;
 306     }
 307     if (ro > 0) {
 308       // Update the Hessian and it's inverse using the BFGS formula
 309 #if 0 // changed by SGJ to compute OBJECTIVE_AND_GRADIENT above
 310       if (axis==-1)
 311         glob.ObjectiveGradient(x_new,g_new,GRADIENT_ONLY);
 312       else {
 313         x_av(axis)=x_new(0);
 314         glob.ObjectiveGradient(x_av,g_av,GRADIENT_ONLY);
 315         g_new(0)=g_av(axis);
 316       }
 317       GC++;
 318       IF_NLOPT_CHECK_EVALS;
 319 #else
 320       if (axis != -1)
 321         g_new(0)=g_av(axis);
 322 #endif
 323
 324       // y=g_new-g
 325       copy(g_new,y);
 326       axpy(-1.0,g,y);
 327
 328       // Check curvature condition
 329       alpha=dot(y,h_dl);
 330       if (alpha <= sqrt(MacEpsilon)*nrm_dl*norm2(y)) {
 331 #ifdef LS_DEBUG
 332         cout << "Curvature condition violated " ;
 333 #endif
 334       }
 335       else {
 336         // Update Hessian
 337         gemv('N',1.0,B,h_dl,0.0,z) ; // z=Bh_dl
 338         beta=-1/dot(h_dl,z) ;
 339         ger(1/alpha,y,y,B) ;
 340         ger(beta,z,z,B) ;
 341
 342         // Update Hessian inverse
 343         gemv('N',1.0,H,y,0.0,z) ; // z=H*y
 344         gemv('T',1.0,H,y,0.0,w) ; // w=y'*H
 345         beta=dot(y,z) ;
 346         beta=(1+beta/alpha)/alpha ;
 347
 348         // It should be possible to do this updating more efficiently, by
 349         // exploiting the fact that (h_dl*y'*H) = transpose(H*y*h_dl')
 350         ger(beta,h_dl,h_dl,H) ;
 351         ger(-1/alpha,z,h_dl,H) ;
 352         ger(-1/alpha,h_dl,w,H) ;
 353       }
 354
 355       if (nrm_dl < norm2(x)*epsilon) {
 356         // Stop criterion (iteration progress) fullfilled
 357 #ifdef LS_DEBUG
 358         cout << "Progress is marginal" ;
 359 #endif
 360         good_enough = 1 ;
 361       }
 362
 363       // Check if we are close to a stationary point located previously
 364       if (box.CloseToMin(x_new, &f_new, eps_cl)) {
 365         // Note that x_new and f_new may be overwritten on exit from CloseToMin
 366 #ifdef LS_DEBUG
 367         cout << "Close to a previously located stationary point, exiting" << endl;
 368 #endif
 369         info = LS_Old ;
 370         good_enough = 1 ;
 371       }
 372
 373       // Update x, g and f
 374       copy(x_new,x) ; copy(g_new,g) ; f=f_new ;
 375
 376 #ifdef LS_DEBUG
 377       cout << " x=" << x << endl ;
 378 #endif
 379
 380     }
 381     else {
 382 #ifdef LS_DEBUG
 383       cout << "Step is no good, ro=" << ro << " delta=" << delta << endl ;
 384 #endif
 385     }
 386
 387   } // wend
 388
 389   // Make sure the routine returns correctly...
 390   // Check if last iterate is outside the boundary
 391   if (box.OutsideBox(x, domain) != 0) {
 392     info=LS_Out; f=DBL_MAX;
 393   }
 394
 395   if (info == LS_Unstable) {
 396     cout << "Local search became unstable. No big deal but exiting anyway\n" ;
 397     exit(1);
 398   }
 399
 400   *mgr=maxgrad ;
 401
 402   T.xvals=x ; T.objval=f ;
 403   if (outside>0)
 404     return LS_Out ;
 405   else
 406     return info ;
 407
 408 #endif
 409 }