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 int local(Trial &T, TBox &box, TBox &domain, double eps_cl, double *mgr,
  23           Global &glob, int axis, RCRVector x_av
  24 #ifdef NLOPT_UTIL_H
  25           , nlopt_stopping *stop
  26 #endif
  27           ) {
  28
  29   int k_max, info, outside ;
  30   int k, i, good_enough, iTmp ;
  31   int n=box.GetDim();
  32
  33   double maxgrad, delta, f, f_new, tmp ;
  34   double alpha, gamma, beta, d2, s2, nom, den, ro ;
  35   double nrm_sd, nrm_hn, snrm_hn, nrm_dl ;
  36   RVector x(n), g(n), h_sd(n), h_dl(n), h_n(n), x_new(n), g_new(n) ;
  37   RVector s(n),y(n),z(n),w(n) ; // Temporary vectors
  38   RMatrix B(n), H(n) ;          // Hessian and it's inverse
  39
  40   k_max = max_iter*n ;
  41   x=T.xvals ;
  42
  43 #ifdef LS_DEBUG
  44   cout << "Local Search, x=" << x << endl;
  45 #endif
  46
  47   if (box.OutsideBox(x, domain) != 0) {
  48     cout << "Starting point is not inside the boundary. Exiting...\n" ;
  49     exit(1) ;
  50     return LS_Out ;
  51   }
  52
  53   // Check if we are close to a stationary point located previously
  54   if (box.CloseToMin(x, &tmp, eps_cl)) {
  55 #ifdef LS_DEBUG
  56      cout << "Close to a previously located stationary point, exiting" << endl;
  57 #endif
  58      T.objval=tmp;
  59      return LS_Old ;
  60    }
  61
  62   // Initially B and H are equal to the identity matrix
  63   B=0 ; H=0 ;
  64   for (i=0 ; i<n ; i++) {
  65     B(i,i)=1 ;
  66     H(i,i)=1 ;
  67   }
  68
  69   RVector g_av(x_av.GetLength());
  70   if (axis==-1) {
  71     f=glob.ObjectiveGradient(x,g,OBJECTIVE_AND_GRADIENT);
  72   }
  73   else {
  74     x_av(axis)=x(0);
  75     f=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
  76     g(0)=g_av(axis);
  77   }
  78   IF_NLOPT_CHECK_EVALS;
  79   FC++;GC++;
  80
  81   if (axis == -1) {
  82     // Skipping AV
  83 #ifdef INI3
  84     // Elaborate scheme to initalize delta
  85     delta=delta_coef*norm2(g) ;
  86     copy(g,z) ;
  87     axpy(1.0,x,z) ;
  88     if (!box.InsideBox(z)) {
  89       if (box.Intersection(x,g,z)==TRUE) {
  90         axpy(-1.0,x,z) ;
  91         delta=min(delta,delta_coef*norm2(z)) ;
  92       }
  93       else {
  94         // Algorithm broke down, use INI1
  95         delta = (1.0/7)*box.ShortestSide(&iTmp) ;
  96       }
  97     }
  98 #endif
  99 #ifdef INI2
 100     // Use INI2 scheme
 101     delta = box.ClosestSide(x)*delta_coef ;
 102     if (delta<MacEpsilon)
 103       // Patch to avoid trust region with radius close to zero
 104       delta = (1.0/7)*box.ShortestSide(&iTmp) ;
 105 #endif
 106 #ifdef INI1
 107     delta = delta_coef*box.ShortestSide(&iTmp) ;
 108 #endif
 109   }
 110   else {
 111     // Use a simple scheme for the 1D minimization (INI1)
 112     delta = (1.0/7.0)*box.ShortestSide(&iTmp) ;
 113   }
 114
 115   k=0 ; good_enough = 0 ; info=LS_New ; outside=0 ;
 116   maxgrad=*mgr ;
 117   while (good_enough == 0) {
 118     k++ ;
 119     if (k>k_max) {
 120 #ifdef LS_DEBUG
 121       cout << "Maximum number of iterations reached\n" ;
 122 #endif
 123       info=LS_MaxIter ;
 124       break ;
 125     }
 126
 127     // Update maximal gradient value
 128     maxgrad=max(maxgrad,normInf(g)) ;
 129
 130     // Steepest descent, h_sd = -g
 131     copy(g,h_sd) ;
 132     scal(-1.0,h_sd) ;
 133     nrm_sd=norm2(h_sd) ;
 134
 135     if (nrm_sd < epsilon) {
 136       // Stop criterion (gradient) fullfilled
 137 #ifdef LS_DEBUG
 138       cout << "Gradient small enough" << endl ;
 139 #endif
 140       good_enough = 1 ;
 141       break ;
 142     }
 143
 144     // Compute Newton step, h_n = -H*g
 145     gemv('N',-1.0, H, g, 0.0, h_n) ;
 146     nrm_hn = norm2(h_n) ;
 147
 148     if (nrm_hn < delta) {
 149       // Pure Newton step
 150       copy(h_n, h_dl) ;
 151 #ifdef LS_DEBUG
 152       cout << "[Newton step]      " ;
 153 #endif
 154     }
 155     else {
 156       gemv('N',1.0,B,g,0.0,z) ;
 157       tmp=dot(g,z) ;
 158       if (tmp==0) {
 159         info = LS_Unstable ;
 160         break ;
 161       }
 162       alpha=(nrm_sd*nrm_sd)/tmp ; // Normalization (N38,eq. 3.30)
 163       scal(alpha,h_sd) ;
 164       nrm_sd=fabs(alpha)*nrm_sd ;
 165
 166       if (nrm_sd >= delta) {
 167         gamma = delta/nrm_sd ; // Normalization (N38, eq. 3.33)
 168         copy(h_sd,h_dl) ;
 169         scal(gamma,h_dl) ;
 170 #ifdef LS_DEBUG
 171         cout << "[Steepest descent]  " ;
 172 #endif
 173       }
 174       else {
 175         // Combination of Newton and SD steps
 176         d2 = delta*delta ;
 177         copy(h_sd,s) ;
 178         s2=nrm_sd*nrm_sd ;
 179         nom = d2 - s2 ;
 180         snrm_hn=nrm_hn*nrm_hn ;
 181         tmp = dot(h_n,s) ;
 182         den = tmp-s2 + sqrt((tmp-d2)*(tmp-d2)+(snrm_hn-d2)*(d2-s2)) ;
 183         if (den==0) {
 184           info = LS_Unstable ;
 185           break ;
 186         }
 187         // Normalization (N38, eq. 3.31)
 188         beta = nom/den ;
 189         copy(h_n,h_dl) ;
 190         scal(beta,h_dl) ;
 191         axpy((1-beta),h_sd,h_dl) ;
 192 #ifdef LS_DEBUG
 193         cout << "[Mixed step]        " ;
 194 #endif
 195       }
 196     }
 197     nrm_dl=norm2(h_dl) ;
 198
 199     //x_new = x+h_dl ;
 200     copy(x,x_new) ;
 201     axpy(1.0,h_dl,x_new) ;
 202
 203     // Check if x_new is inside the box
 204     iTmp=box.OutsideBox(x_new, domain) ;
 205     if (iTmp == 1) {
 206 #ifdef LS_DEBUG
 207       cout << "x_new is outside the box " << endl ;
 208 #endif
 209       outside++ ;
 210       if (outside>max_outside_steps) {
 211         // Previous point was also outside, exit
 212         break ;
 213       }
 214     }
 215     else if (iTmp == 2) {
 216 #ifdef LS_DEBUG
 217       cout << " x_new is outside the domain" << endl ;
 218 #endif
 219       info=LS_Out ;
 220       break ;
 221     }
 222     else {
 223       outside=0 ;
 224     }
 225
 226     // Compute the gain
 227     if (axis==-1)
 228       f_new=glob.ObjectiveGradient(x_new,g_new,OBJECTIVE_AND_GRADIENT);
 229     else {
 230       x_av(axis)=x_new(0);
 231       f_new=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
 232     }
 233     IF_NLOPT_CHECK_EVALS;
 234     FC++; GC++;
 235     gemv('N',0.5,B,h_dl,0.0,z);
 236     ro = (f_new-f) / (dot(g,h_dl) + dot(h_dl,z)); // Quadratic model
 237     if (ro > 0.75) {
 238       delta = delta*2;
 239     }
 240     if (ro < 0.25) {
 241       delta = delta/3;
 242     }
 243     if (ro > 0) {
 244       // Update the Hessian and it's inverse using the BFGS formula
 245 #if 0 // changed by SGJ to compute OBJECTIVE_AND_GRADIENT above
 246       if (axis==-1)
 247         glob.ObjectiveGradient(x_new,g_new,GRADIENT_ONLY);
 248       else {
 249         x_av(axis)=x_new(0);
 250         glob.ObjectiveGradient(x_av,g_av,GRADIENT_ONLY);
 251         g_new(0)=g_av(axis);
 252       }
 253       GC++;
 254       IF_NLOPT_CHECK_EVALS;
 255 #else
 256       if (axis != -1)
 257         g_new(0)=g_av(axis);
 258 #endif
 259
 260       // y=g_new-g
 261       copy(g_new,y);
 262       axpy(-1.0,g,y);
 263
 264       // Check curvature condition
 265       alpha=dot(y,h_dl);
 266       if (alpha <= sqrt(MacEpsilon)*nrm_dl*norm2(y)) {
 267 #ifdef LS_DEBUG
 268         cout << "Curvature condition violated " ;
 269 #endif
 270       }
 271       else {
 272         // Update Hessian
 273         gemv('N',1.0,B,h_dl,0.0,z) ; // z=Bh_dl
 274         beta=-1/dot(h_dl,z) ;
 275         ger(1/alpha,y,y,B) ;
 276         ger(beta,z,z,B) ;
 277
 278         // Update Hessian inverse
 279         gemv('N',1.0,H,y,0.0,z) ; // z=H*y
 280         gemv('T',1.0,H,y,0.0,w) ; // w=y'*H
 281         beta=dot(y,z) ;
 282         beta=(1+beta/alpha)/alpha ;
 283
 284         // It should be possible to do this updating more efficiently, by
 285         // exploiting the fact that (h_dl*y'*H) = transpose(H*y*h_dl')
 286         ger(beta,h_dl,h_dl,H) ;
 287         ger(-1/alpha,z,h_dl,H) ;
 288         ger(-1/alpha,h_dl,w,H) ;
 289       }
 290
 291       if (nrm_dl < norm2(x)*epsilon) {
 292         // Stop criterion (iteration progress) fullfilled
 293 #ifdef LS_DEBUG
 294         cout << "Progress is marginal" ;
 295 #endif
 296         good_enough = 1 ;
 297       }
 298
 299       // Check if we are close to a stationary point located previously
 300       if (box.CloseToMin(x_new, &f_new, eps_cl)) {
 301         // Note that x_new and f_new may be overwritten on exit from CloseToMin
 302 #ifdef LS_DEBUG
 303         cout << "Close to a previously located stationary point, exiting" << endl;
 304 #endif
 305         info = LS_Old ;
 306         good_enough = 1 ;
 307       }
 308
 309       // Update x, g and f
 310       copy(x_new,x) ; copy(g_new,g) ; f=f_new ;
 311
 312 #ifdef LS_DEBUG
 313       cout << " x=" << x << endl ;
 314 #endif
 315
 316     }
 317     else {
 318 #ifdef LS_DEBUG
 319       cout << "Step is no good, ro=" << ro << " delta=" << delta << endl ;
 320 #endif
 321     }
 322
 323   } // wend
 324
 325   // Make sure the routine returns correctly...
 326   // Check if last iterate is outside the boundary
 327   if (box.OutsideBox(x, domain) != 0) {
 328     info=LS_Out; f=DBL_MAX;
 329   }
 330
 331   if (info == LS_Unstable) {
 332     cout << "Local search became unstable. No big deal but exiting anyway\n" ;
 333     exit(1);
 334   }
 335
 336   T.xvals=x ; T.objval=f ;
 337   *mgr=maxgrad ;
 338   if (outside>0)
 339     return LS_Out ;
 340   else
 341     return info ;
 342 }